Nucleation Experiment

Author

Holl et al.

Published

August 27, 2025

show code
knitr::opts_chunk$set(echo = TRUE)

After first submission

  • Data and R codes available: github, then it will be submited to zenodo

  • Inclusion of species diversity of recruiting plants: item I.C, with rarefaction and Proportion of Biotic/Abiotic Seed Dispersal and Planted/Unplanted species.

  • Table with statistical details from the models (file "itatinga_NUC_ModelSummaries.xlsx")

  • Fixed LiDAR figure (SE of the mean)

  • Fixed plots with empty recruits

  • Considering only small recruits in species richness rarefaction curves

  • Analyses separating Baccharis + Vernonanthura

Setup environment

1) R Environment

show code
path <- getwd()
version
               _                                
platform       x86_64-w64-mingw32               
arch           x86_64                           
os             mingw32                          
crt            ucrt                             
system         x86_64, mingw32                  
status                                          
major          4                                
minor          3.3                              
year           2024                             
month          02                               
day            29                               
svn rev        86002                            
language       R                                
version.string R version 4.3.3 (2024-02-29 ucrt)
nickname       Angel Food Cake                  
show code
if(!require(pacman)){install.packages("pacman")}
Carregando pacotes exigidos: pacman
show code
pkg_list <-c("tidyverse", "lme4", "DHARMa", "emmeans", "multcomp", "iNEXT",
             "car", "performance", "cowplot", "MuMIn", "ggpattern")
pacman::p_load(pkg_list, character.only = TRUE)

2) Read data

trees_clean: data of trees (planted/not-planted)

show code
trees_clean <- readxl::read_xlsx(paste0(path, "/Itatinga_Nuc_trees.xlsx"))
trees_clean <- trees_clean |>
  dplyr::mutate(scientificName = dplyr::if_else(is.na(scientificName), "Unclassified", scientificName))

grass_cover: grass cover per plot (%)

show code
grass_cover <- readxl::read_xlsx(paste0(path, "/Itatinga_Nuc_grasscover.xlsx"))

lidar: Canopy Height Metric (CHM_mean, CHM_SD, and CHM_CV), Canopy cover (Canopy_cover), Rugosity (Rugosity_SD, and Rugosity_CV), Leaf Area Index (LAI_mean, LAI_SD, and LAI_CV), and LAIunder (LAIunder_mean, LAIunder_SD, LAIunder_CV).

show code
lidar2022_wide <- readxl::read_xlsx(
  paste0(path, "/tab_metrics_LiDAR_2022_v2.xlsx")) |>
  dplyr::select(Trat:LAIunder_SD) |>
  tidyr::separate(Trat, into = c("Bloco", "Perc_cobertura","Tratamento"), sep = "_") |>
  dplyr::mutate(Tratamento = dplyr::if_else(is.na(Tratamento), "controle", Tratamento),
                Tratamento = paste(Perc_cobertura, Tratamento, sep="_"),
                ano =  2022) |>
  dplyr::select(-Perc_cobertura) 
New names:
• `` -> `...1`
Warning: Expected 3 pieces. Missing pieces filled with `NA` in 5 rows [5, 10, 15, 20,
25].
show code
lidar2024_wide <- readxl::read_xlsx(
  paste0(path, "/tab_metrics_LiDAR_2024_v2.xlsx"))|>
  dplyr::select(Trat:LAIunder_SD) |>
  tidyr::separate(Trat, into = c("Bloco", "Perc_cobertura","Tratamento"), sep = "_") |>
  dplyr::mutate(Tratamento = dplyr::if_else(is.na(Tratamento), "controle", Tratamento),
                Tratamento = paste(Perc_cobertura, Tratamento, sep="_"),
                ano =  2024) |>
  dplyr::select(-Perc_cobertura)
New names:
• `` -> `...1`
Warning: Expected 3 pieces. Missing pieces filled with `NA` in 5 rows [5, 10, 15, 20,
25].
show code
lidar_wide <- dplyr::bind_rows(lidar2022_wide, lidar2024_wide) |> # 40 x 12
   dplyr::mutate(
     Tratamento = factor(Tratamento, 
                         levels = c("25_tabuleiro", "25_faixa",  "50_tabuleiro", "50_faixa", "100_controle"),
                         labels = c("Nuclei 25%", "Strips 25%", "Nuclei 50%", "Strips 50%", "Planting 100%")))

Table - List of sampled species (planted and recruits)

two tables:

show code
trees_clean |> dplyr::select(RP, family, scientificName) |>
  unique()|>
  dplyr::mutate(valor=TRUE)|>
  tidyr::pivot_wider(names_from = RP, values_from = valor, values_fill = FALSE) |>
  dplyr::arrange(desc(planted), family, scientificName) # |> writexl::write_xlsx("tab1.xlsx")
# A tibble: 108 × 4
   family        scientificName                      planted regenerant
   <chr>         <chr>                               <lgl>   <lgl>     
 1 Anacardiaceae Astronium graveolens                TRUE    FALSE     
 2 Annonaceae    Annona cacans                       TRUE    TRUE      
 3 Bignoniaceae  Handroanthus impetiginosus          TRUE    FALSE     
 4 Bignoniaceae  Handroanthus ochraceus              TRUE    TRUE      
 5 Boraginaceae  Cordia sellowiana                   TRUE    TRUE      
 6 Euphorbiaceae Croton floribundus                  TRUE    TRUE      
 7 Euphorbiaceae Croton urucurana                    TRUE    FALSE     
 8 Euphorbiaceae Mabea fistulifera                   TRUE    TRUE      
 9 Fabaceae      Anadenanthera colubrina  var. cebil TRUE    TRUE      
10 Fabaceae      Apuleia leiocarpa                   TRUE    FALSE     
# ℹ 98 more rows
show code
trees_clean |> dplyr::select(RP, family, scientificName) |>
  unique() |> with(table(RP))
RP
   planted regenerant 
        37         86 

I. Results - Field data

A) Survival of planted trees

(i) Organizing data

show code
trees_clean |> with(table(RP, scientificName)) |>
  as.data.frame() |> 
  dplyr::filter(Freq>0) |> # writexl::write_xlsx("tab_spp_plantio.xlsx")
  dplyr::filter(RP=="planted")
        RP                      scientificName Freq
1  planted Anadenanthera colubrina  var. cebil  204
2  planted                       Annona cacans    4
3  planted                   Apuleia leiocarpa   11
4  planted                Astronium graveolens    5
5  planted              Cariniana estrellensis    2
6  planted               Cecropia pachystachya  117
7  planted                    Cedrela fissilis   28
8  planted                      Ceiba speciosa   54
9  planted             Citharexylum myrianthum   16
10 planted              Copaifera langsdorffii   24
11 planted                   Cordia sellowiana   10
12 planted                  Croton floribundus  831
13 planted                    Croton urucurana  726
14 planted       Enterolobium contortisiliquum   77
15 planted                   Erythrina mulungu   32
16 planted                    Eugenia uniflora    9
17 planted                   Ficus guaranitica   82
18 planted               Gallesia integrifolia   50
19 planted                    Genipa americana    3
20 planted                   Guazuma ulmifolia  323
21 planted          Handroanthus impetiginosus   30
22 planted              Handroanthus ochraceus   20
23 planted                  Hymenaea courbaril   12
24 planted                         Inga edulis   85
25 planted                           Inga vera  741
26 planted                    Lafoensia pacari   34
27 planted                   Luehea divaricata   48
28 planted                   Mabea fistulifera   20
29 planted              Machaerium acutifolium   17
30 planted                Myroxylon peruiferum    2
31 planted                  Peltophorum dubium   29
32 planted                 Platypodium elegans    7
33 planted                   Psidium myrtoides    2
34 planted                    Pterogyne nitens    3
35 planted                Senegalia polyphylla  474
36 planted                   Senna macranthera    1
37 planted                        Unclassified  372

Reading scientific names

show code
survival_planted2024 <- trees_clean |> # 5416 stems
  dplyr::filter(ano==2024)|> 
  dplyr::select(bloco:arv, Hm, scientificName, RP) |> # 3,199 stems
  unique() |> # 1487
  dplyr::filter(RP=="planted") |> # 1135
  dplyr::group_by(tratamento,bloco, parcela)|>
  dplyr::summarise(n = dplyr::n(), .groups = "drop")|>
  dplyr::arrange(desc(n))

(ii) Survival analysis

24 seedlings per plot

show code
survival_planted2024 |>
  with(table(bloco, parcela, by=tratamento))
, , by = 100

     parcela
bloco 1 2 3 4
    1 1 1 1 1
    2 1 1 1 1
    3 1 1 1 1
    4 1 1 1 1
    5 1 1 1 1

, , by = F25

     parcela
bloco 1 2 3 4
    1 1 1 1 1
    2 1 1 1 1
    3 1 1 1 1
    4 1 1 1 1
    5 1 1 1 1

, , by = F50

     parcela
bloco 1 2 3 4
    1 1 1 1 1
    2 1 1 1 1
    3 1 1 1 1
    4 1 1 1 1
    5 1 1 1 1

, , by = N25

     parcela
bloco 1 2 3 4
    1 1 1 1 1
    2 1 1 1 1
    3 1 1 1 1
    4 1 1 1 1
    5 1 1 1 1

, , by = N50

     parcela
bloco 1 2 3 4
    1 1 1 1 1
    2 1 1 1 1
    3 1 1 1 1
    4 1 1 1 1
    5 1 1 1 1
show code
survival_planted2024 |> dplyr::arrange(tratamento, bloco, parcela)
# A tibble: 100 × 4
   tratamento bloco parcela     n
   <chr>      <chr> <chr>   <int>
 1 100        1     1          17
 2 100        1     2          15
 3 100        1     3          12
 4 100        1     4          19
 5 100        2     1          11
 6 100        2     2          11
 7 100        2     3          17
 8 100        2     4          13
 9 100        3     1          13
10 100        3     2          16
# ℹ 90 more rows
show code
mod_surv1 <- lme4::glmer(cbind(n, 24 - n) ~ tratamento + (1|bloco), data = survival_planted2024, family = binomial(link = "logit"))
mod_surv_null <- lme4::glmer(cbind(n, 24 - n) ~ 0 +(1|bloco), data=survival_planted2024, family = binomial(link = "logit"))
AIC(mod_surv1, mod_surv_null)
              df      AIC
mod_surv1      6 554.2907
mod_surv_null  1 564.5393
show code
par(mfrow=c(1,2))
DHARMa::plotResiduals(mod_surv1)
DHARMa::plotQQunif(mod_surv1)

Marginal \(R^2\) / Conditional \(R^2\)

show code
MuMIn::r.squaredGLMM(mod_surv1)
Warning: the null model is only correct if all the variables it uses are identical 
to those used in fitting the original model.
                  R2m       R2c
theoretical 0.1291593 0.4869551
delta       0.1160836 0.4376569
show code
0.129/ 0.487
[1] 0.2648871

Which it is 0.265

show code
sjPlot::tab_model(mod_surv1)
  cbind(n,24-n)
Predictors Odds Ratios CI p
(Intercept) 1.10 0.79 – 1.52 0.580
tratamento [F25] 0.80 0.62 – 1.03 0.090
tratamento [F50] 0.91 0.70 – 1.18 0.472
tratamento [N25] 0.58 0.45 – 0.75 <0.001
tratamento [N50] 0.84 0.65 – 1.09 0.192
Random Effects
σ2 3.29
τ00 bloco 0.10
ICC 0.03
N bloco 5
Observations 100
Marginal R2 / Conditional R2 0.010 / 0.038
show code
broom.mixed::tidy(mod_surv1) |> dplyr::mutate(dplyr::across(4:7, round, 3)) 
Registered S3 method overwritten by 'broom':
  method        from 
  nobs.multinom MuMIn
Warning: There was 1 warning in `dplyr::mutate()`.
ℹ In argument: `dplyr::across(4:7, round, 3)`.
Caused by warning:
! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
Supply arguments directly to `.fns` through an anonymous function instead.

  # Previously
  across(a:b, mean, na.rm = TRUE)

  # Now
  across(a:b, \(x) mean(x, na.rm = TRUE))
# A tibble: 6 × 7
  effect   group term            estimate std.error statistic p.value
  <chr>    <chr> <chr>              <dbl>     <dbl>     <dbl>   <dbl>
1 fixed    <NA>  (Intercept)        0.092     0.166     0.553   0.58 
2 fixed    <NA>  tratamentoF25     -0.222     0.131    -1.70    0.09 
3 fixed    <NA>  tratamentoF50     -0.094     0.131    -0.719   0.472
4 fixed    <NA>  tratamentoN25     -0.544     0.132    -4.12    0    
5 fixed    <NA>  tratamentoN50     -0.171     0.131    -1.31    0.192
6 ran_pars bloco sd__(Intercept)    0.309    NA        NA      NA    
show code
summary(mod_surv1)
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: cbind(n, 24 - n) ~ tratamento + (1 | bloco)
   Data: survival_planted2024

      AIC       BIC    logLik -2*log(L)  df.resid 
    554.3     569.9    -271.1     542.3        94 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.0363 -0.8992 -0.1404  1.0640  3.6310 

Random effects:
 Groups Name        Variance Std.Dev.
 bloco  (Intercept) 0.0956   0.3092  
Number of obs: 100, groups:  bloco, 5

Fixed effects:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)    0.09194    0.16634   0.553   0.5805    
tratamentoF25 -0.22197    0.13080  -1.697   0.0897 .  
tratamentoF50 -0.09389    0.13068  -0.719   0.4724    
tratamentoN25 -0.54411    0.13224  -4.115 3.88e-05 ***
tratamentoN50 -0.17069    0.13072  -1.306   0.1916    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) trtF25 trtF50 trtN25
tratamntF25 -0.393                     
tratamntF50 -0.393  0.500              
tratamntN25 -0.389  0.494  0.495       
tratamntN50 -0.393  0.500  0.500  0.495
show code
exp(-0.544)
[1] 0.5804219
show code
(52.3/(100-52.3))/ (39.2/(100-39.2))
[1] 1.700595

Planted trees in treatment F25 have 19.91% less likely to survive (OR=exp(-0.222)) than in treatment 100% (\(p<0.1\))

Planted trees in treatment N25 have 41.96% less likely to survive (OR=exp(-0.544)) than in treatment 100% (\(p<0.001\))

show code
round((1 - exp(-0.544))*100,2)
[1] 41.96
show code
emmeans::emmeans(mod_surv1, ~ tratamento, type = "response") |> 
  multcomp::cld(details = TRUE, adjust="sidak")
$emmeans
 tratamento  prob     SE  df asymp.LCL asymp.UCL .group
 N25        0.389 0.0398 Inf     0.293     0.495  1    
 F25        0.468 0.0414 Inf     0.364     0.574  12   
 N50        0.480 0.0415 Inf     0.376     0.586   2   
 F50        0.500 0.0416 Inf     0.394     0.605   2   
 100        0.523 0.0415 Inf     0.417     0.627   2   

Confidence level used: 0.95 
Conf-level adjustment: sidak method for 5 estimates 
Intervals are back-transformed from the logit scale 
P value adjustment: sidak method for 10 tests 
Tests are performed on the log odds ratio scale 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

$comparisons
 contrast  odds.ratio    SE  df null z.ratio p.value
 F25 / N25       1.38 0.183 Inf    1   2.436  0.1390
 N50 / N25       1.45 0.192 Inf    1   2.825  0.0463
 N50 / F25       1.05 0.138 Inf    1   0.392  1.0000
 F50 / N25       1.57 0.207 Inf    1   3.407  0.0066
 F50 / F25       1.14 0.149 Inf    1   0.980  0.9810
 F50 / N50       1.08 0.141 Inf    1   0.588  0.9997
 100 / N25       1.72 0.228 Inf    1   4.115  0.0004
 100 / F25       1.25 0.163 Inf    1   1.697  0.6093
 100 / N50       1.19 0.155 Inf    1   1.306  0.8809
 100 / F50       1.10 0.144 Inf    1   0.719  0.9983

P value adjustment: sidak method for 10 tests 
Tests are performed on the log odds ratio scale 

(iii) Table S2 - Survival

show code
survival_planted2024 |>
  dplyr::group_by( tratamento)|>
  dplyr::summarise(Mean = round(mean((n/24)*100),1),
                  Max = round(max((n/24)*100),1),
                  Min = round(min((n/24)*100),1) ) |>
  dplyr::mutate(values = paste0(Mean," (",Min," - ",Max,")")) |>
  dplyr::select(-c(Mean, Max, Min)) 
# A tibble: 5 × 2
  tratamento values            
  <chr>      <chr>             
1 100        52.3 (29.2 - 79.2)
2 F25        46.9 (12.5 - 70.8)
3 F50        50 (29.2 - 75)    
4 N25        39.2 (8.3 - 79.2) 
5 N50        48.1 (20.8 - 75)  

B) Grass cover

(i) Organizing data

How planting scheme affected grass cover inside and outside planted areas?

Excluding treatment 100:

show code
ordem <- c("N25", "F25","N50", "F50", "100")
grass_cover <- readxl::read_xlsx(
  paste0(path, "/Itatinga_Nuc_grasscover.xlsx"))|>
  dplyr::mutate(
    grass_perc = dplyr::if_else(gramineas_perc>=1, 0.9999, gramineas_perc),
    tratamento = ordered(tratamento, 
                         levels = ordem,
                         labels = c("Nuclei 25%", "Strips 25%", "Nuclei 50%", "Strips 50%", "Planting 100%"))) |>
   droplevels()
grass_cover2024 <- grass_cover |> dplyr::filter(tratamento!="Planting 100%") |> dplyr::filter(ano==2024)
grass_cover2024 |>
  with(table(bloco, parcela, by=tratamento))
, , by = Nuclei 25%

     parcela
bloco 1 2 3 4
    1 2 2 2 2
    2 2 2 2 2
    3 2 2 2 2
    4 2 2 2 2
    5 2 2 2 2

, , by = Strips 25%

     parcela
bloco 1 2 3 4
    1 2 2 2 2
    2 2 2 2 2
    3 2 2 2 2
    4 2 2 2 2
    5 2 2 2 2

, , by = Nuclei 50%

     parcela
bloco 1 2 3 4
    1 2 2 2 2
    2 2 2 2 2
    3 2 2 2 2
    4 2 2 2 2
    5 2 2 2 2

, , by = Strips 50%

     parcela
bloco 1 2 3 4
    1 2 2 2 2
    2 2 2 2 2
    3 2 2 2 2
    4 2 2 2 2
    5 2 2 2 2

, , by = Planting 100%

     parcela
bloco 1 2 3 4
    1 0 0 0 0
    2 0 0 0 0
    3 0 0 0 0
    4 0 0 0 0
    5 0 0 0 0

Too many zeros? Zero-Inflated model? There is few zeros. Only binomial fit

show code
grass_cover2024 |>
  ggplot(aes(x=grass_perc))+
  geom_histogram()+
  facet_grid(plantio~tratamento)+theme_classic()
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

(ii) Analysis

show code
grass_cover2024 |>
  with(table(bloco, parcela, by=tratamento))
, , by = Nuclei 25%

     parcela
bloco 1 2 3 4
    1 2 2 2 2
    2 2 2 2 2
    3 2 2 2 2
    4 2 2 2 2
    5 2 2 2 2

, , by = Strips 25%

     parcela
bloco 1 2 3 4
    1 2 2 2 2
    2 2 2 2 2
    3 2 2 2 2
    4 2 2 2 2
    5 2 2 2 2

, , by = Nuclei 50%

     parcela
bloco 1 2 3 4
    1 2 2 2 2
    2 2 2 2 2
    3 2 2 2 2
    4 2 2 2 2
    5 2 2 2 2

, , by = Strips 50%

     parcela
bloco 1 2 3 4
    1 2 2 2 2
    2 2 2 2 2
    3 2 2 2 2
    4 2 2 2 2
    5 2 2 2 2

, , by = Planting 100%

     parcela
bloco 1 2 3 4
    1 0 0 0 0
    2 0 0 0 0
    3 0 0 0 0
    4 0 0 0 0
    5 0 0 0 0
show code
grass_cover2024 <-grass_cover2024 |>
  dplyr::mutate(Tratamento=factor(tratamento, ordered = FALSE))
modelo_gc24a <- glmmTMB::glmmTMB(grass_perc ~ Tratamento*plantio + (1 | bloco),
                        family = glmmTMB::beta_family(), data = grass_cover2024)
modelo_gc24b <- grass_cover2024 |> 
  dplyr::mutate(tratamento=factor(tratamento,ordered = FALSE))|>
  with(glmmTMB::glmmTMB(grass_perc ~ tratamento+plantio + (1 | bloco),
                   ziformula = ~ 1, family = glmmTMB::beta_family()))
Warning in glmmTMB::glmmTMB(grass_perc ~ tratamento + plantio + (1 | bloco), :
use of the 'data' argument is recommended
show code
modelo_gc24c <-  grass_cover2024 |> 
  dplyr::mutate(tratamento=factor(tratamento,ordered = FALSE))|>
  with(glmmTMB::glmmTMB(grass_perc ~ plantio + (1 | bloco), 
                                  family = glmmTMB::beta_family()))
Warning in glmmTMB::glmmTMB(grass_perc ~ plantio + (1 | bloco), family =
glmmTMB::beta_family()): use of the 'data' argument is recommended
show code
modelo_gc24d <-  grass_cover2024 |> 
  dplyr::mutate(tratamento=factor(tratamento,ordered = FALSE))|>
  with(glmmTMB::glmmTMB(grass_perc ~ tratamento + (1 | bloco), 
                                  family = glmmTMB::beta_family()))
Warning in glmmTMB::glmmTMB(grass_perc ~ tratamento + (1 | bloco), family =
glmmTMB::beta_family()): use of the 'data' argument is recommended
show code
modelo_gc24null <-  grass_cover2024 |> 
  dplyr::mutate(tratamento=factor(tratamento,ordered = FALSE))|>
  with(glmmTMB::glmmTMB(grass_perc ~ 0 + (1 | bloco), 
                        family = glmmTMB::beta_family()))
Warning in glmmTMB::glmmTMB(grass_perc ~ 0 + (1 | bloco), family =
glmmTMB::beta_family()): use of the 'data' argument is recommended
show code
AIC(modelo_gc24a, modelo_gc24b, modelo_gc24c, modelo_gc24d, modelo_gc24null)
                df       AIC
modelo_gc24a    10 -61.25388
modelo_gc24b     8 -53.84400
modelo_gc24c     4 -51.58877
modelo_gc24d     6 -21.66566
modelo_gc24null  2 -10.55579

Analysis of residuals

show code
par(mfrow=c(1,2))
DHARMa::plotQQunif(modelo_gc24a)
DHARMa::plotResiduals(modelo_gc24a)

Is there any difference between treatments considering inside/outside?

show code
resultado_24 <- emmeans::emmeans(modelo_gc24a, ~ Tratamento|plantio) |> 
  multcomp::cld( details = TRUE)
resultado_24$emmeans
plantio = CP:
 Tratamento emmean    SE  df asymp.LCL asymp.UCL .group
 Strips 50% -0.727 0.228 Inf    -1.173    -0.280  1    
 Strips 25% -0.300 0.224 Inf    -0.739     0.138  12   
 Nuclei 25%  0.138 0.223 Inf    -0.299     0.575   23  
 Nuclei 50%  0.708 0.228 Inf     0.262     1.154    3  

plantio = SP:
 Tratamento emmean    SE  df asymp.LCL asymp.UCL .group
 Strips 50%  0.937 0.231 Inf     0.485     1.389  1    
 Nuclei 50%  0.945 0.231 Inf     0.493     1.397  1    
 Nuclei 25%  1.036 0.232 Inf     0.581     1.490  1    
 Strips 25%  1.166 0.234 Inf     0.708     1.625  1    

Results are given on the logit (not the response) scale. 
Confidence level used: 0.95 
Results are given on the log odds ratio (not the response) scale. 
P value adjustment: tukey method for comparing a family of 4 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Answer: yes.

Is there any difference between inside/outside considering treatments?

show code
resultado_24a <- emmeans::emmeans(modelo_gc24a, ~ plantio|Tratamento) |> 
  multcomp::cld(details = TRUE)
resultado_24a$emmeans
Tratamento = Nuclei 25%:
 plantio emmean    SE  df asymp.LCL asymp.UCL .group
 CP       0.138 0.223 Inf    -0.299     0.575  1    
 SP       1.036 0.232 Inf     0.581     1.490   2   

Tratamento = Strips 25%:
 plantio emmean    SE  df asymp.LCL asymp.UCL .group
 CP      -0.300 0.224 Inf    -0.739     0.138  1    
 SP       1.166 0.234 Inf     0.708     1.625   2   

Tratamento = Nuclei 50%:
 plantio emmean    SE  df asymp.LCL asymp.UCL .group
 CP       0.708 0.228 Inf     0.262     1.154  1    
 SP       0.945 0.231 Inf     0.493     1.397  1    

Tratamento = Strips 50%:
 plantio emmean    SE  df asymp.LCL asymp.UCL .group
 CP      -0.727 0.228 Inf    -1.173    -0.280  1    
 SP       0.937 0.231 Inf     0.485     1.389   2   

Results are given on the logit (not the response) scale. 
Confidence level used: 0.95 
Results are given on the log odds ratio (not the response) scale. 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Answer: yes.

show code
broom.mixed::tidy(modelo_gc24a) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))|> knitr::kable()
effect component group term estimate std.error statistic p.value
fixed cond NA (Intercept) 0.138 0.223 0.620 0.535
fixed cond NA TratamentoStrips 25% -0.439 0.316 -1.388 0.165
fixed cond NA TratamentoNuclei 50% 0.569 0.318 1.790 0.073
fixed cond NA TratamentoStrips 50% -0.865 0.319 -2.710 0.007
fixed cond NA plantioSP 0.897 0.321 2.794 0.005
fixed cond NA TratamentoStrips 25%:plantioSP 0.569 0.453 1.256 0.209
fixed cond NA TratamentoNuclei 50%:plantioSP -0.660 0.454 -1.454 0.146
fixed cond NA TratamentoStrips 50%:plantioSP 0.766 0.454 1.688 0.091
ran_pars cond bloco sd__(Intercept) 0.000 NA NA NA
show code
sjPlot::tab_model(modelo_gc24a)
  grass perc
Predictors Estimates CI p
(Intercept) 1.15 0.74 – 1.78 0.535
Tratamento [Strips 25%] 0.64 0.35 – 1.20 0.165
Tratamento [Nuclei 50%] 1.77 0.95 – 3.30 0.073
Tratamento [Strips 50%] 0.42 0.23 – 0.79 0.007
plantio [SP] 2.45 1.31 – 4.60 0.005
Tratamento [Strips 25%] ×
plantio [SP]
1.77 0.73 – 4.30 0.209
Tratamento [Nuclei 50%] ×
plantio [SP]
0.52 0.21 – 1.26 0.146
Tratamento [Strips 50%] ×
plantio [SP]
2.15 0.88 – 5.24 0.091
Random Effects
σ2 0.16
τ00 bloco 0.00
N bloco 5
Observations 160
Marginal R2 / Conditional R2 0.736 / NA
show code
lme4::VarCorr(modelo_gc24a)  # ou summary(modelo_gc24a)$varcor

Conditional model:
 Groups Name        Std.Dev.  
 bloco  (Intercept) 1.0464e-05
show code
MuMIn::r.squaredGLMM(modelo_gc24a)
Warning: Can't compute random effect variances. Some variance components equal
  zero. Your model may suffer from singularity (see `?lme4::isSingular`
  and `?performance::check_singularity`).
  Decrease the `tolerance` level to force the calculation of random effect
  variances, or impose priors on your random effects parameters (using
  packages like `brms` or `glmmTMB`).
           R2m       R2c
[1,] 0.7355586 0.7355586
show code
m_fix <- update(modelo_gc24a, . ~ . - (1|bloco))
anova(modelo_gc24a, m_fix) 
Data: grass_cover2024
Models:
m_fix: grass_perc ~ Tratamento + plantio + Tratamento:plantio, zi=~0, disp=~1
modelo_gc24a: grass_perc ~ Tratamento * plantio + (1 | bloco), zi=~0, disp=~1
             Df     AIC     BIC logLik deviance Chisq Chi Df Pr(>Chisq)
m_fix         9 -63.254 -35.577 40.627  -81.254                        
modelo_gc24a 10 -61.254 -30.502 40.627  -81.254     0      1          1
show code
performance::check_singularity(modelo_gc24a)
[1] TRUE

(iii) Figure - Grass cover

Small letters: Comparison between treatments, only inside planted areas.

Capital letters: Comparison between treatments, only outside planted areas.

show code
letras <- data.frame(
  tratamento = rep(levels(as.factor(grass_cover$tratamento)),each=2),
  plantio = levels(as.factor(grass_cover$plantio)),
  label = c("BC      ", "      a",
            "AB      ", "      a",
            "C      ", "      a",
            "A      ", "      a",
            "A",""),  
  y = c(0.6, 0.82, 
        0.5, 0.8, 
        0.65, 0.8, 
        0.35, 0.8,
        0.35, .1)
  ) |>
  dplyr::mutate(
    tratamento = forcats::fct_recode(tratamento,
                                     "Nuclei\n25%" = "Nuclei 25%",
                                     "Strips\n25%" = "Strips 25%",
                                     "Nuclei\n50%" = "Nuclei 50%",
                                     "Strips\n50%" = "Strips 50%",
                                     "Planting\n100%" = "Planting 100%"
                                     )
  )

fig_grass <- grass_cover |>  
  dplyr::filter(ano==2024) |>
  dplyr::group_by(tratamento, plantio) |>
  dplyr::summarise(grass_cover_mean = mean(gramineas_perc),
                   sd = sd(gramineas_perc),
                   n = dplyr::n(),
                   se = sd/sqrt(n),
                   .groups="drop") |>
  dplyr::mutate(
    tratamento = forcats::fct_recode(tratamento,
                                     "Nuclei\n25%" = "Nuclei 25%",
                                     "Strips\n25%" = "Strips 25%",
                                     "Nuclei\n50%" = "Nuclei 50%",
                                     "Strips\n50%" = "Strips 50%",
                                     "Planting\n100%" = "Planting 100%"
                                     )
  ) |>
  ggplot(aes(x = tratamento, y = grass_cover_mean, fill = plantio)) +
  geom_bar(position = position_dodge2(preserve = "single"), 
           stat = "identity", color = "black") +
  geom_errorbar(aes(ymin = grass_cover_mean - se, ymax = grass_cover_mean + se),
                width = .2, position = position_dodge(0.9)) +
  scale_fill_manual(
    labels = c('Planted', 'Unplanted'), 
    values = c("#92D050", "white")) +
  geom_text(data = letras, aes(x = tratamento, y = y, label = label), size=4, vjust = -0.5) +
  labs(fill = "", 
       y = "Grass cover (%)",
       x = "") +
  scale_y_continuous(expand = c(0,0), 
                     labels = scales::percent_format(accuracy = 1), 
                     limits = c(0,1)) +
  theme_classic()+
  theme(
    axis.text.x = element_text(color="black", hjust = 0.5, size=11),
    axis.text.y = element_text(color="black", size=11),
    axis.title.y = element_text(color="black", size=13),
    legend.position = c(0.84, 0.99),
    legend.text = element_text(size = 11, margin = margin(0,0,0, l=1, "mm")),
    legend.key.size = unit(3, "mm"),
    legend.spacing = unit(0, "mm"),
    legend.background = element_blank())
Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.
show code
fig_grass

show code
# ggsave("fig_grass.png", plot=fig_grass, width=90, height = 90 ,unit="mm")

(iv) Table S2 - Grass cover

show code
grass_cover |>
  dplyr::group_by(ano, tratamento, plantio)|>
  dplyr::summarise(Mean = mean(gramineas_perc)*100,
                  Max = max(gramineas_perc)*100,
                  Min = min(gramineas_perc)*100) |>
  dplyr::mutate(values = paste0(Mean," (",Min," - ",Max,")")) |>
  dplyr::select(-c(Mean, Max, Min)) |>
  tidyr::pivot_wider(names_from = tratamento, values_from = values)
`summarise()` has grouped output by 'ano', 'tratamento'. You can override using
the `.groups` argument.
# A tibble: 4 × 7
# Groups:   ano [2]
    ano plantio `Nuclei 25%`       `Strips 25%`      `Nuclei 50%`   `Strips 50%`
  <dbl> <chr>   <chr>              <chr>             <chr>          <chr>       
1  2022 CP      29.2 (0 - 85)      16.9 (0 - 100)    34.9 (0 - 100) 12.35 (0 - …
2  2022 SP      65.25 (20 - 100)   65.175 (23 - 100) 51.6 (0 - 100) 60.35 (0 - …
3  2024 CP      54.07 (8.2 - 95)   39.91 (4 - 95)    56.3485 (2.64… 29.9795 (5 …
4  2024 SP      76.835 (39 - 97.3) 75.45 (29 - 100)  73.627 (26 - … 77.9725 (49…
# ℹ 1 more variable: `Planting 100%` <chr>

C) Recruitment

Only recruits (smaller and larger), not including planted individuals

show code
species_table <- readxl::read_xlsx("itatinga_Nuc_regen_species_abundance.xlsx")
species_table|>
  dplyr::select(Treatment, Block, Plot, Planted)|>
  unique()|>
  with(table(Block, Plot, by=Treatment))
, , by = Full_Planted

     Plot
Block P1 P2 P3 P4
   B1  1  1  1  1
   B2  1  1  1  1
   B3  1  1  1  1
   B4  1  1  1  1
   B5  1  1  1  1

, , by = Nuclei_25

     Plot
Block P1 P2 P3 P4
   B1  2  2  2  2
   B2  2  2  2  2
   B3  2  2  2  2
   B4  2  2  2  2
   B5  2  2  2  2

, , by = Nuclei_50

     Plot
Block P1 P2 P3 P4
   B1  2  2  2  2
   B2  2  2  2  2
   B3  2  2  2  2
   B4  2  2  2  2
   B5  2  2  2  2

, , by = Strip_25

     Plot
Block P1 P2 P3 P4
   B1  2  2  2  2
   B2  2  2  2  2
   B3  2  2  2  2
   B4  2  2  2  2
   B5  2  2  2  2

, , by = Strip_50

     Plot
Block P1 P2 P3 P4
   B1  2  2  2  2
   B2  2  2  2  2
   B3  2  2  2  2
   B4  2  2  2  2
   B5  2  2  2  2

Viewing top species

show code
total_smaller <- sum(species_table$n_sapling)
total_larger <- sum(species_table$n_tree)
species_table |>
  dplyr::group_by(scientificName)|>
  dplyr::summarise(
    sum_large = sum(n_tree, na.rm=T)/total_larger,
    sum_small = sum(n_sapling, na.rm=T)/total_smaller,
  )|>
  dplyr::arrange(desc(sum_small)) |>
  head()
# A tibble: 6 × 3
  scientificName             sum_large sum_small
  <chr>                          <dbl>     <dbl>
1 Baccharis dracunculifolia    0.0740     0.489 
2 Vernonanthura polyanthes     0          0.313 
3 Solanum sisymbriifolium      0.00157    0.0388
4 Croton floribundus           0.00157    0.0384
5 Senegalia polyphylla         0.00472    0.0305
6 Solanum granuloso-leprosum   0.0205     0.0233

(i) Considering Baccharis + Vernonanthura

Only small regenerants

show code
dados_list <- species_table |>
  dplyr::group_by(Treatment, Planted, scientificName) |>
  dplyr::summarise(n_total = sum(n_sapling), .groups = "drop") |>
  dplyr::group_by(Treatment, Planted) |>
  dplyr::summarise(abundancias = list(n_total), .groups = "drop") |>
  dplyr::mutate(nome_grupo = paste(Treatment, Planted, sep = "_")) %>% 
  { setNames(.$abundancias, .$nome_grupo) }
res_inext <- iNEXT::iNEXT(
  dados_list,
  q = 0,
  size=seq(0, 2000, 100),
  datatype = "abundance"
)

estimativas <- res_inext$iNextEst |>
  dplyr::bind_rows() |>
  tidyr::separate(Assemblage, into = c("Treatment", "Percent", "Planted"), sep = "_")

Figure species richness

Figure Recruits’ species richness

For Hill q = 0

Only small regenerants and considering ruderal species (Vernonanthura & Baccharis)

show code
dados_list <- species_table |>
  dplyr::group_by(Treatment, Planted, scientificName) |>
  dplyr::summarise(n_total = sum(n_sapling), .groups = "drop") |>
  dplyr::group_by(Treatment, Planted) |>
  dplyr::summarise(abundancias = list(n_total), .groups = "drop") |>
  dplyr::mutate(nome_grupo = paste(Treatment, Planted, sep = "_")) %>% 
  { setNames(.$abundancias, .$nome_grupo) }
res_inext <- iNEXT::iNEXT(
  dados_list,
  q = 0,
  size=seq(0, 2000, 100),
  datatype = "abundance"
)

estimativas <- res_inext$iNextEst |>
  dplyr::bind_rows() |>
  tidyr::separate(Assemblage, into = c("Treatment", "Percent", "Planted"), sep = "_")

curvas_df0 <- estimativas |> 
  dplyr::filter(Order.q==0)|> 
  dplyr::filter(!is.na(SC.LCL))|>
  dplyr::mutate(
    Percent = ifelse(Percent=="Planted", "100", Percent),
    Percent = ordered(Percent, levels = c("25", "50", "100")),
    Planted = ifelse(Planted=="Inside", "Planted", "Unplanted"),
    Treatment = ifelse(Treatment=="Full", "Planting 100%", Treatment),
    Treatment = factor(Treatment, levels = c("Nuclei", "Strip", "Planting 100%"), ordered = TRUE)
  ) 
fig_recruit_rich <- curvas_df0 |>
  ggplot(aes(x = m, y = qD, color = Treatment, linetype = Planted)) +
  geom_line(linewidth = 1) +
  geom_ribbon(aes(ymin = qD.LCL, ymax = qD.UCL, fill = Treatment), alpha = 0.2, colour = NA) +
   geom_point(
    data = dplyr::filter(curvas_df0, Method == "Observed"),
    size = 2.5
  ) +
  facet_grid(Percent~., 
             labeller =  as_labeller(
               c("25" = "(a) Planting 25%",
                 "50" = "(b) Planting 50%",
                 "100" = "(c) Planting 100%")
               )) +
  scale_color_manual(values=c("red", "blue", "darkgreen"))+
  scale_fill_manual(values=c("red", "blue", "darkgreen"))+
  labs(
    x = "Number of sampling units (regenerants)",
    y = "Species diversity (Hill Index, q=0)"
    ) +
  ggthemes::theme_clean()+
  scale_x_continuous(expand = c(0.001, 0), breaks = seq(0, 3500, 500)) +
  scale_y_continuous(expand = c(0, 0))+
  theme(
     panel.grid.major.y = element_blank(),
    panel.grid.major.x = element_blank(),
    axis.text.x = element_text(color = "black", hjust = 1, size = 10),
    axis.text.y = element_text(color = "black", size = 9),
    axis.title.y = element_text(color = "black", size = 11),
    strip.text = element_text(color = "black", hjust = 0.5, size = 10),
    strip.placement = "outside",
    legend.text = element_text(size = 11),
    legend.title = element_blank(),
    legend.key.size = unit(5, "mm"),
    legend.box.background = element_blank(),
    legend.box.margin = margin(t = 5, r = 0, b = 0, l = 0, "mm"),
    legend.background = element_blank(),
    legend.margin = margin(0, 0, 0, 0),
    legend.position = c(0.25, 0.25),
    panel.border = element_rect(colour = "black", fill = NA),
    panel.spacing = unit(0, "lines")
  )
fig_recruit_rich

show code
# ggsave("fig_richness_all_small_recruits.png", plot=fig_recruit_rich, width=90, height = 210 ,unit="mm")

Table S2 - Recruitment (all)

Estimating Density of recruits (\(ind.ha^{-1}\)) and considering ruderal species (Baccharis and Vernonanthura)

show code
density_recruits_n_m2 <- species_table |> 
  dplyr::mutate(
    n_small = dplyr::if_else(is.na(n_sapling), 0, n_sapling),
    n_large = dplyr::if_else(is.na(n_tree), 0, n_tree)) |>
  dplyr::group_by(Treatment, Block, Planted) |>
  dplyr::summarise(
    n_small_m2 = (sum(n_small, na.rm=T))/(144*4), 
    n_large_m2 = (sum(n_large, na.rm=T))/(144*4), 
    .groups = "drop")

dens_ha <- density_recruits_n_m2 |>
  dplyr::mutate(
    n_small_ha = n_small_m2 * 10000,
    n_large_ha = n_large_m2 * 10000
  )
dens_P_UP <- dens_ha |>
  tidyr::pivot_wider(
    names_from = Planted,
    values_from = c(n_small_ha, n_large_ha),
    names_sep = "_"
  )

prop_table <- tidyr::tibble(
  Treatment = unique(sort(species_table$Treatment)),
  prop_P = c(1, 0.25, 0.5, 0.25, 0.5),
  prop_UP = 1 - prop_P
)

result_regenerant_all <- dens_P_UP |>
  dplyr::left_join(prop_table, by = "Treatment") |>
  dplyr::mutate(
    n_small_total_ha = dplyr::coalesce(n_small_ha_Inside * prop_P, 0)  + dplyr::coalesce(n_small_ha_Outside * prop_UP, 0),
    n_large_total_ha = dplyr::coalesce(n_large_ha_Inside * prop_P, 0) + dplyr::coalesce(n_large_ha_Outside * prop_UP, 0)
  ) |>
  dplyr::select(Treatment,
         starts_with("n_small_ha"),
         starts_with("n_large_ha"),
         n_small_total_ha,
         n_large_total_ha
  )

Smaller recruits

  • Only in planted areas
show code
result_regenerant_all |>
  dplyr::group_by(Treatment)|>
  dplyr::summarise(
    mean = round(mean(n_small_ha_Inside, na.rm=TRUE),1),
    min = round(min(n_small_ha_Inside, na.rm=TRUE),0),
    max = round(max(n_small_ha_Inside, na.rm=TRUE),0)
  )|>
  dplyr::mutate(
    larger = sprintf("%.1f (%.1f - %.1f)", 
                                mean, 
                                min, 
                                max)) |>
  dplyr::select(Treatment, larger) |>
  tidyr::pivot_wider(
    names_from = Treatment,
    values_from = larger
    )
# A tibble: 1 × 5
  Full_Planted             Nuclei_25                Nuclei_50  Strip_25 Strip_50
  <chr>                    <chr>                    <chr>      <chr>    <chr>   
1 3048.6 (2274.0 - 4705.0) 4871.5 (3490.0 - 7101.0) 6107.6 (2… 4816.0 … 6024.3 …
  • Only in unplanted areas
show code
result_regenerant_all |>
  dplyr::group_by(Treatment)|>
  dplyr::summarise(
    mean = round(mean(n_small_ha_Outside, na.rm=TRUE), 1),
    min = round(min(n_small_ha_Outside, na.rm=TRUE), 0),
    max = round(max(n_small_ha_Outside, na.rm=TRUE), 0)
  )|>
  dplyr::mutate(
    larger = sprintf("%.1f (%.1f - %.1f)", 
                                mean, 
                                min, 
                                max)) |>
  dplyr::select(Treatment, larger) |>
  tidyr::pivot_wider(
    names_from = Treatment,
    values_from = larger
    )
Warning: There were 2 warnings in `dplyr::summarise()`.
The first warning was:
ℹ In argument: `min = round(min(n_small_ha_Outside, na.rm = TRUE), 0)`.
ℹ In group 1: `Treatment = "Full_Planted"`.
Caused by warning in `min()`:
! nenhum argumento não faltante para min; retornando Inf
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
# A tibble: 1 × 5
  Full_Planted     Nuclei_25                Nuclei_50          Strip_25 Strip_50
  <chr>            <chr>                    <chr>              <chr>    <chr>   
1 NaN (Inf - -Inf) 3961.8 (2344.0 - 8247.0) 3461.8 (2031.0 - … 3392.4 … 3871.5 …
  • In total (planted + unplanted areas)
show code
result_regenerant_all |>
  dplyr::group_by(Treatment)|>
  dplyr::summarise(
    mean = round(mean(n_small_total_ha, na.rm=TRUE),1),
    min = round(min(n_small_total_ha, na.rm=TRUE),0),
    max = round(max(n_small_total_ha, na.rm=TRUE),0)
  )|>
  dplyr::mutate(
    larger = sprintf("%.1f (%.1f - %.1f)", 
                                mean, 
                                min, 
                                max)) |>
  dplyr::select(Treatment, larger) |>
  tidyr::pivot_wider(
    names_from = Treatment,
    values_from = larger
    )
# A tibble: 1 × 5
  Full_Planted             Nuclei_25               Nuclei_50   Strip_25 Strip_50
  <chr>                    <chr>                   <chr>       <chr>    <chr>   
1 3048.6 (2274.0 - 4705.0) 2094.6 (872.0 - 6185.0) 2392.4 (10… 1874.1 … 2474.0 …

Larger recruits

All treatments and inside/outside areas showed similar values

  • Only in planted areas
show code
result_regenerant_all |>
  dplyr::group_by(Treatment)|>
  dplyr::summarise(
    mean = round(mean(n_large_ha_Inside, na.rm=TRUE),1),
    min = round(min(n_large_ha_Inside, na.rm=TRUE),0),
    max = round(max(n_large_ha_Inside, na.rm=TRUE),0)
  )|>
  dplyr::mutate(
    larger = sprintf("%.1f (%.1f - %.1f)", 
                                mean, 
                                min, 
                                max)) |>
  dplyr::select(Treatment, larger) |>
  tidyr::pivot_wider(
    names_from = Treatment,
    values_from = larger
    )
# A tibble: 1 × 5
  Full_Planted          Nuclei_25           Nuclei_50          Strip_25 Strip_50
  <chr>                 <chr>               <chr>              <chr>    <chr>   
1 222.2 (104.0 - 486.0) 229.2 (0.0 - 608.0) 163.2 (35.0 - 295… 152.8 (… 225.7 (…
  • Only in unplanted areas
show code
result_regenerant_all |>
  dplyr::group_by(Treatment)|>
  dplyr::summarise(
    mean = round(mean(n_large_ha_Outside, na.rm=TRUE),1),
    min = round(min(n_large_ha_Outside, na.rm=TRUE),0),
    max = round(max(n_large_ha_Outside, na.rm=TRUE),0)
  )|>
  dplyr::mutate(
    larger = sprintf("%.1f (%.1f - %.1f)", 
                                mean, 
                                min, 
                                max)) |>
  dplyr::select(Treatment, larger) |>
  tidyr::pivot_wider(
    names_from = Treatment,
    values_from = larger
    )
Warning: There were 2 warnings in `dplyr::summarise()`.
The first warning was:
ℹ In argument: `min = round(min(n_large_ha_Outside, na.rm = TRUE), 0)`.
ℹ In group 1: `Treatment = "Full_Planted"`.
Caused by warning in `min()`:
! nenhum argumento não faltante para min; retornando Inf
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
# A tibble: 1 × 5
  Full_Planted     Nuclei_25            Nuclei_50             Strip_25  Strip_50
  <chr>            <chr>                <chr>                 <chr>     <chr>   
1 NaN (Inf - -Inf) 364.6 (17.0 - 851.0) 163.2 (122.0 - 208.0) 302.1 (1… 381.9 (…
  • In total (planted + unplanted areas)
show code
result_regenerant_all |>
  dplyr::group_by(Treatment)|>
  dplyr::summarise(
    mean = round(mean(n_large_total_ha, na.rm = TRUE), 1),
    min = round(min(n_large_total_ha, na.rm = TRUE), 0),
    max = round(max(n_large_total_ha, na.rm = TRUE), 0)
  )|>
  dplyr::mutate(
    larger = sprintf("%.1f (%.1f - %.1f)", 
                                mean, 
                                min, 
                                max)) |>
  dplyr::select(Treatment, larger) |>
  tidyr::pivot_wider(
    names_from = Treatment,
    values_from = larger
    )
# A tibble: 1 × 5
  Full_Planted          Nuclei_25           Nuclei_50          Strip_25 Strip_50
  <chr>                 <chr>               <chr>              <chr>    <chr>   
1 222.2 (104.0 - 486.0) 165.4 (0.0 - 638.0) 81.6 (17.0 - 148.… 132.4 (… 151.9 (…

Analysis for smaller recruits abundance

Data:

show code
small_df <- species_table |>
  dplyr::group_by(Treatment, Block, Plot, Planted)|>
  dplyr::summarise(n = sum(n_sapling), .groups="drop")

small_df |>
  ggplot(aes(x=Planted, y=n, color=Block))+
  geom_point()+facet_wrap(~Treatment)+theme_bw()

Models:

show code
modelo_reg1 <- glmmTMB::glmmTMB(n ~ Treatment*Planted + (1 | Block),
                   # family = poisson(),
                   family = glmmTMB::nbinom2(),
                   data = small_df)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
show code
modelo_reg2 <-  glmmTMB::glmmTMB(n ~ Treatment + (1 | Block),
                   # family = poisson(),
                   family = glmmTMB::nbinom2(),
                   data = small_df)
modelo_reg3 <-  glmmTMB::glmmTMB(n ~ Planted + (1 | Block),
                   # family = poisson(),
                   family = glmmTMB::nbinom2(),
                   data = small_df)
modelo_reg_null <-  glmmTMB::glmmTMB(n ~ 0 + (1 | Block),
                   # family = poisson(), # Poisson
                   family = glmmTMB::nbinom2(), # NB
                   data = small_df)

# NBinom 1782.688 (null) |  1734.985 (Treatment*Planted) | 1756.665 (only Treament) | 1743.097  (only Planted)
# Poisson  4032.177 |   3373.854        |   3813.769        |   3759.763

AIC(modelo_reg_null, modelo_reg1, modelo_reg2, modelo_reg3)
                df      AIC
modelo_reg_null  2 1782.688
modelo_reg1     11 1734.985
modelo_reg2      7 1756.665
modelo_reg3      4 1743.097

non-positive-definite Hessian matrix indicates that there are problems in estimating the model parameters. This may be related to redundancies in the random effects, perfect separation, few data for some combinations of factors, or problems in choosing the distribution. Therefore, let’s check the singularity:

show code
performance::check_singularity(modelo_reg1)
[1] FALSE
show code
performance::check_model(modelo_reg1)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
`check_outliers()` does not yet support models of class `glmmTMB`.

Exploring residuals

show code
par(mfrow=c(2,2))
DHARMa::plotResiduals(modelo_reg1)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
show code
DHARMa::plotQQunif(modelo_reg1)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
show code
DHARMa::plotResiduals(modelo_reg3)
DHARMa::plotQQunif(modelo_reg3)

show code
car::Anova(modelo_reg1)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: n
                    Chisq Df Pr(>Chisq)    
Treatment         23.4883  4  0.0001011 ***
Planted           30.8344  1   2.81e-08 ***
Treatment:Planted  1.8521  3  0.6036618    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Is there any difference between treatments considering only inside or outside?

show code
emmeans::emmeans(modelo_reg1, ~ Treatment|Planted) |> 
  multcomp::cld(adjust="bonferroni")
Planted = Inside:
 Treatment    emmean    SE  df asymp.LCL asymp.UCL .group
 Full_Planted   3.76 0.175 Inf      3.31      4.21  1    
 Strip_25       4.24 0.173 Inf      3.79      4.68   2   
 Nuclei_25      4.24 0.174 Inf      3.80      4.69   2   
 Nuclei_50      4.40 0.173 Inf      3.96      4.85   2   
 Strip_50       4.42 0.173 Inf      3.98      4.87   2   

Planted = Outside:
 Treatment    emmean    SE  df asymp.LCL asymp.UCL .group
 Strip_25       3.82 0.175 Inf      3.38      4.25  1    
 Nuclei_50      3.82 0.175 Inf      3.39      4.26  1    
 Nuclei_25      3.96 0.174 Inf      3.52      4.39  1    
 Strip_50       3.98 0.174 Inf      3.55      4.42  1    
 Full_Planted nonEst    NA  NA        NA        NA       

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for varying numbers of estimates 
P value adjustment: bonferroni method for varying numbers of tests 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Answer: Yes (inside) and no (outside).

Is there any difference between treatments not considering inside/outside?

show code
emmeans::emmeans(modelo_reg2, ~ Treatment) |> 
    multcomp::cld(adjust="bonferroni")
 Treatment    emmean    SE  df asymp.LCL asymp.UCL .group
 Full_Planted   3.76 0.177 Inf      3.30      4.22  1    
 Strip_25       4.05 0.155 Inf      3.65      4.45  12   
 Nuclei_25      4.11 0.155 Inf      3.71      4.51  12   
 Nuclei_50      4.16 0.155 Inf      3.76      4.56  12   
 Strip_50       4.23 0.154 Inf      3.83      4.63   2   

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 5 estimates 
P value adjustment: bonferroni method for 10 tests 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Answer: No.

Is there any difference between inside/outside considering treatments?

show code
emmeans::emmeans(modelo_reg1, ~ Planted|Treatment) |> 
    multcomp::cld(adjust="bonferroni")
Treatment = Full_Planted:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Inside    3.76 0.175 Inf      3.42      4.10  1    
 Outside nonEst    NA  NA        NA        NA       

Treatment = Nuclei_25:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Outside   3.96 0.174 Inf      3.57      4.35  1    
 Inside    4.24 0.174 Inf      3.85      4.63  1    

Treatment = Nuclei_50:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Outside   3.82 0.175 Inf      3.43      4.21  1    
 Inside    4.40 0.173 Inf      4.02      4.79   2   

Treatment = Strip_25:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Outside   3.82 0.175 Inf      3.42      4.21  1    
 Inside    4.24 0.173 Inf      3.85      4.63   2   

Treatment = Strip_50:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Outside   3.98 0.174 Inf      3.59      4.37  1    
 Inside    4.42 0.173 Inf      4.03      4.81   2   

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for varying numbers of estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Answer: Yes.

Table S3

show code
sjPlot::tab_model(modelo_reg1)
Model matrix is rank deficient. Parameters
  `TreatmentStrip_50:PlantedOutside` were not estimable.
  n
Predictors Incidence Rate Ratios CI p
(Intercept) 42.99 30.52 – 60.55 <0.001
Treatment [Nuclei_25] 1.62 1.19 – 2.20 0.002
Treatment [Nuclei_50] 1.90 1.40 – 2.59 <0.001
Treatment [Strip_25] 1.61 1.19 – 2.19 0.002
Treatment [Strip_50] 1.94 1.43 – 2.63 <0.001
Planted [Outside] 0.64 0.48 – 0.87 0.005
Treatment [Nuclei_25] ×
Planted [Outside]
1.17 0.76 – 1.80 0.478
Treatment [Nuclei_50] ×
Planted [Outside]
0.87 0.56 – 1.33 0.515
Treatment [Strip_25] ×
Planted [Outside]
1.02 0.66 – 1.57 0.937
Random Effects
σ2 0.21
τ00 Block 0.09
ICC 0.30
N Block 5
Observations 180
Marginal R2 / Conditional R2 0.164 / 0.412
show code
broom.mixed::tidy(modelo_reg1) |>
  dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3))) |> dplyr::select(-c(group,component))
# A tibble: 11 × 6
   effect   term                            estimate std.error statistic p.value
   <chr>    <chr>                              <dbl>     <dbl>     <dbl>   <dbl>
 1 fixed    (Intercept)                        3.76      0.175    21.5     0    
 2 fixed    TreatmentNuclei_25                 0.481     0.156     3.08    0.002
 3 fixed    TreatmentNuclei_50                 0.644     0.156     4.13    0    
 4 fixed    TreatmentStrip_25                  0.476     0.156     3.05    0.002
 5 fixed    TreatmentStrip_50                  0.662     0.156     4.25    0    
 6 fixed    PlantedOutside                    -0.439     0.155    -2.83    0.005
 7 fixed    TreatmentNuclei_25:PlantedOuts…    0.156     0.219     0.709   0.478
 8 fixed    TreatmentNuclei_50:PlantedOuts…   -0.143     0.22     -0.652   0.515
 9 fixed    TreatmentStrip_25:PlantedOutsi…    0.017     0.22      0.079   0.937
10 fixed    TreatmentStrip_50:PlantedOutsi…   NA        NA        NA      NA    
11 ran_pars sd__(Intercept)                    0.301    NA        NA      NA    

Analysis for larger recruits abundance

Data:

show code
large_regenerants_2024 <- species_table |>
  dplyr::group_by(Treatment, Block, Plot, Planted)|>
  dplyr::summarise(n = sum(n_tree), .groups="drop")

Plotting

show code
large_regenerants_2024 |>
  ggplot(aes(x=Planted, y=n, color=Block))+
  geom_point()+facet_wrap(~Treatment)+theme_bw()

Models

show code
modelo_reg1 <- glmmTMB::glmmTMB(n ~ Treatment*Planted + (1 | Block),
                   ziformula = ~1,
                    # family = poisson(),
                     family = glmmTMB::nbinom2(),
                    # family = glmmTMB::truncated_nbinom2(),
                   data = large_regenerants_2024)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
show code
modelo_reg2 <-  glmmTMB::glmmTMB(n ~ Treatment + (1 | Block),
                   ziformula = ~1,
                   #  family = poisson(),
                    family = glmmTMB::nbinom2(),
                  #  family = glmmTMB::truncated_nbinom2(),
                   data = large_regenerants_2024)
modelo_reg3 <-  glmmTMB::glmmTMB(n ~ Planted + (1 | Block),
                   ziformula = ~1,
                   # family = poisson(),
                    family = glmmTMB::nbinom2(),
                   #  family = glmmTMB::truncated_nbinom2(),
                   data = large_regenerants_2024)
modelo_reg_null <-  glmmTMB::glmmTMB(n ~ 0 + (1 | Block),
                   ziformula = ~1,
                    # family = poisson(), # ZIP
                     family = glmmTMB::nbinom2(), # ZINB
                    # family = glmmTMB::truncated_nbinom2(), # hurdle
                   data = large_regenerants_2024)
# Model (null) | (treat*planted) | (treat) | (planted)
# ZINB 885.72       |   869.63       |  868.93  |   864.21
# ZIP  1058.14   | 1015.06       |  1025.38 |   1024.49
# Hurdle 880.75 | 872.21        |   868.74      |   975.88

AIC(modelo_reg_null, modelo_reg1, modelo_reg2, modelo_reg3)
                df      AIC
modelo_reg_null  3 885.7221
modelo_reg1     12 869.6264
modelo_reg2      8 868.9342
modelo_reg3      5 864.2072
show code
performance::check_model(modelo_reg1)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
`check_outliers()` does not yet support models of class `glmmTMB`.

Exploring residuals

show code
par(mfrow=c(2,2))
DHARMa::plotResiduals(modelo_reg1)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
show code
DHARMa::plotQQunif(modelo_reg1)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
show code
DHARMa::plotResiduals(modelo_reg3)
DHARMa::plotQQunif(modelo_reg3)

show code
MuMIn::r.squaredGLMM(modelo_reg3)
Warning: the null model is only correct if all the variables it uses are identical 
to those used in fitting the original model.
                 R2m        R2c
delta     0.03380077 0.03380078
lognormal 0.05135140 0.05135140
trigamma  0.01854201 0.01854201
show code
# 0.03253385/0.03665597

0.0328/ 0.0338 variance explained by fixed effects and by the entire model (around 0.9999997 % is explained by fixed effects)

show code
performance::r2(modelo_reg1)
Warning: Can't compute random effect variances. Some variance components equal
  zero. Your model may suffer from singularity (see `?lme4::isSingular`
  and `?performance::check_singularity`).
  Decrease the `tolerance` level to force the calculation of random effect
  variances, or impose priors on your random effects parameters (using
  packages like `brms` or `glmmTMB`).
[1] NA
show code
broom.mixed::tidy(modelo_reg1) |>
  dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3))) |>
  dplyr::select(-c(component, group))|>
  knitr::kable()
effect term estimate std.error statistic p.value
fixed (Intercept) 1.175 0.251 4.690 0.000
fixed TreatmentNuclei_25 0.064 0.369 0.173 0.863
fixed TreatmentNuclei_50 -0.286 0.369 -0.776 0.438
fixed TreatmentStrip_25 -0.353 0.369 -0.957 0.339
fixed TreatmentStrip_50 0.028 0.354 0.080 0.936
fixed PlantedOutside 0.510 0.346 1.476 0.140
fixed TreatmentNuclei_25:PlantedOutside -0.045 0.492 -0.092 0.927
fixed TreatmentNuclei_50:PlantedOutside -0.515 0.507 -1.018 0.309
fixed TreatmentStrip_25:PlantedOutside 0.169 0.499 0.339 0.734
fixed TreatmentStrip_50:PlantedOutside NA NA NA NA
fixed (Intercept) -3.510 3.051 -1.150 0.250
ran_pars sd__(Intercept) 0.000 NA NA NA
show code
car::Anova(modelo_reg1)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: n
                   Chisq Df Pr(>Chisq)  
Treatment         7.3216  4    0.11984  
Planted           5.5546  1    0.01843 *
Treatment:Planted 1.9111  3    0.59106  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Is there any difference between treatments considering inside/outside?

show code
emmeans::emmeans(modelo_reg1, ~ Treatment|Planted) |> 
  multcomp::cld(adjust = "bonferroni")
Planted = Inside:
 Treatment    emmean    SE  df asymp.LCL asymp.UCL .group
 Strip_25      0.822 0.284 Inf    0.0899      1.55  1    
 Nuclei_50     0.889 0.284 Inf    0.1589      1.62  1    
 Full_Planted  1.175 0.251 Inf    0.5298      1.82  1    
 Strip_50      1.203 0.260 Inf    0.5327      1.87  1    
 Nuclei_25     1.239 0.287 Inf    0.4988      1.98  1    

Planted = Outside:
 Treatment    emmean    SE  df asymp.LCL asymp.UCL .group
 Nuclei_50     0.884 0.277 Inf    0.1920      1.58  1    
 Strip_25      1.501 0.259 Inf    0.8532      2.15  1    
 Nuclei_25     1.704 0.276 Inf    1.0145      2.39  1    
 Strip_50      1.714 0.235 Inf    1.1262      2.30  1    
 Full_Planted nonEst    NA  NA        NA        NA       

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for varying numbers of estimates 
P value adjustment: bonferroni method for varying numbers of tests 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Answer: No.

Is there any difference between inside/outside considering treatments?

show code
emmeans::emmeans(modelo_reg1, ~ Planted|Treatment) |> 
  multcomp::cld(adjust = "bonferroni")
Treatment = Full_Planted:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Inside   1.175 0.251 Inf     0.684      1.67  1    
 Outside nonEst    NA  NA        NA        NA       

Treatment = Nuclei_25:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Inside   1.239 0.287 Inf     0.595      1.88  1    
 Outside  1.704 0.276 Inf     1.085      2.32  1    

Treatment = Nuclei_50:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Outside  0.884 0.277 Inf     0.263      1.51  1    
 Inside   0.889 0.284 Inf     0.254      1.52  1    

Treatment = Strip_25:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Inside   0.822 0.284 Inf     0.185      1.46  1    
 Outside  1.501 0.259 Inf     0.920      2.08  1    

Treatment = Strip_50:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Inside   1.203 0.260 Inf     0.620      1.79  1    
 Outside  1.714 0.235 Inf     1.186      2.24  1    

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for varying numbers of estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Answer: No.

Propotion of smaller recruits from planted species

  • With ruderal species

Getting proportions of planted/unplanted species

show code
prop_regenerant_full <- species_table |>
  dplyr::mutate(
    planted_tree = ifelse(`Planting group`=="NA", "No", "Yes"))|>
  dplyr::filter(n_sapling>0)|>
  dplyr::group_by(Treatment, Planted, Block, Plot) |>
  dplyr::summarise(
    n_species = dplyr::n(),
    n_ind = sum(n_sapling),
    n_planted_ind = sum(ifelse(planted_tree == "Yes", n_sapling, 0), na.rm = TRUE),
    n_planted_ind = ifelse(is.na(n_planted_ind), 0, n_planted_ind),
    n_planted_spp = sum(planted_tree == "Yes"), 
    n_planted_spp = ifelse(is.na(n_planted_spp), 0, n_planted_spp),
    .groups = "drop") |>
  dplyr::mutate(
    prop_planted_ind = n_planted_ind / n_ind,
    prop_planted_spp = n_planted_spp / n_species
  )

prop_regenerant_full |>
  dplyr::group_by(Treatment, Planted)|>
  dplyr::summarise(
    prop_planted_mean = mean(prop_planted_ind, na.rm=T),
    prop_planted_min = min(prop_planted_ind, na.rm=T),
    prop_planted_max = max(prop_planted_ind, na.rm=T),
    prop_planted_spp_mean = mean(prop_planted_spp, na.rm=T),
    prop_planted_spp_min = min(prop_planted_spp, na.rm=T),
    prop_planted_spp_max = max(prop_planted_spp, na.rm=T),
    .groups = "drop"
  ) |>
  dplyr::mutate(
    prop_planted = sprintf("%.1f (%.1f - %.1f)", 
                                round(prop_planted_mean,3)*100, 
                                round(prop_planted_min,3)*100, 
                                round(prop_planted_max,3)*100),
    prop_planted_spp = sprintf("%.1f (%.1f - %.1f)", 
                                    round(prop_planted_spp_mean,3)*100, 
                                    round(prop_planted_spp_min,3)*100, 
                                    round(prop_planted_spp_max,3)*100)
    ) |>
  dplyr::select(Treatment,Planted, prop_planted, prop_planted_spp) |>
  tidyr::pivot_wider(
    names_from = Treatment,
    values_from = c(prop_planted, prop_planted_spp)
    ) # |> writexl::write_xlsx("resulta_prop_recruit.xlsx")
# A tibble: 2 × 11
  Planted prop_planted_Full_Plan…¹ prop_planted_Nuclei_25 prop_planted_Nuclei_50
  <chr>   <chr>                    <chr>                  <chr>                 
1 Inside  10.5 (0.0 - 43.3)        3.9 (0.0 - 14.8)       7.7 (0.0 - 44.4)      
2 Outside <NA>                     10.4 (0.0 - 91.2)      10.0 (0.0 - 69.2)     
# ℹ abbreviated name: ¹​prop_planted_Full_Planted
# ℹ 7 more variables: prop_planted_Strip_25 <chr>, prop_planted_Strip_50 <chr>,
#   prop_planted_spp_Full_Planted <chr>, prop_planted_spp_Nuclei_25 <chr>,
#   prop_planted_spp_Nuclei_50 <chr>, prop_planted_spp_Strip_25 <chr>,
#   prop_planted_spp_Strip_50 <chr>

Statistical analysis

show code
mod_binomial_full <-lme4::glmer(
  cbind(n_planted_ind, n_ind - n_planted_ind) ~ Treatment*Planted + (1 | Block),
  family = binomial(link = "logit"),
  data = prop_regenerant_full)
fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
show code
performance::check_overdispersion(mod_binomial_full) # no
# Overdispersion test

 dispersion ratio = 1.748
          p-value = 0.272
No overdispersion detected.
show code
performance::check_model(mod_binomial_full)

show code
DHARMa::plotResiduals(mod_binomial_full)

show code
emmeans::emmeans(mod_binomial_full, pairwise ~ Treatment | Planted) |> multcomp::cld(adjust="sidak")
Planted = Inside:
 Treatment    emmean    SE  df asymp.LCL asymp.UCL .group
 Nuclei_25     -3.35 0.292 Inf     -4.10    -2.598  1    
 Strip_50      -2.63 0.271 Inf     -3.32    -1.933   2   
 Nuclei_50     -2.59 0.269 Inf     -3.28    -1.901   2   
 Strip_25      -2.40 0.269 Inf     -3.09    -1.708   2   
 Full_Planted  -2.30 0.279 Inf     -3.02    -1.583   2   

Planted = Outside:
 Treatment    emmean    SE  df asymp.LCL asymp.UCL .group
 Strip_25      -2.79 0.285 Inf     -3.50    -2.084  1    
 Nuclei_25     -2.48 0.278 Inf     -3.17    -1.789  12   
 Nuclei_50     -2.30 0.276 Inf     -2.99    -1.610   2   
 Strip_50      -1.56 0.265 Inf     -2.22    -0.903    3  
 Full_Planted nonEst    NA  NA        NA        NA       

Results are given on the logit (not the response) scale. 
Confidence level used: 0.95 
Conf-level adjustment: sidak method for varying numbers of estimates 
Results are given on the log odds ratio (not the response) scale. 
P value adjustment: sidak method for varying numbers of tests 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
show code
emmeans::emmeans(mod_binomial_full, pairwise ~  Planted | Treatment) |> multcomp::cld(adjust="sidak")
Treatment = Full_Planted:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Inside   -2.30 0.279 Inf     -2.85    -1.753  1    
 Outside nonEst    NA  NA        NA        NA       

Treatment = Nuclei_25:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Inside   -3.35 0.292 Inf     -4.00    -2.695  1    
 Outside  -2.48 0.278 Inf     -3.10    -1.860   2   

Treatment = Nuclei_50:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Inside   -2.59 0.269 Inf     -3.19    -1.991  1    
 Outside  -2.30 0.276 Inf     -2.92    -1.681   2   

Treatment = Strip_25:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Outside  -2.79 0.285 Inf     -3.43    -2.156  1    
 Inside   -2.40 0.269 Inf     -3.00    -1.797   2   

Treatment = Strip_50:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Inside   -2.63 0.271 Inf     -3.23    -2.023  1    
 Outside  -1.56 0.265 Inf     -2.16    -0.971   2   

Results are given on the logit (not the response) scale. 
Confidence level used: 0.95 
Conf-level adjustment: sidak method for varying numbers of estimates 
Results are given on the log odds ratio (not the response) scale. 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
show code
broom.mixed::tidy(mod_binomial_full) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))
# A tibble: 10 × 7
   effect   group term                      estimate std.error statistic p.value
   <chr>    <chr> <chr>                        <dbl>     <dbl>     <dbl>   <dbl>
 1 fixed    <NA>  (Intercept)                 -2.30      0.279    -8.24    0    
 2 fixed    <NA>  TreatmentNuclei_25          -1.05      0.188    -5.56    0    
 3 fixed    <NA>  TreatmentNuclei_50          -0.293     0.151    -1.95    0.051
 4 fixed    <NA>  TreatmentStrip_25           -0.101     0.151    -0.665   0.506
 5 fixed    <NA>  TreatmentStrip_50           -0.329     0.154    -2.14    0.033
 6 fixed    <NA>  PlantedOutside               1.06      0.128     8.33    0    
 7 fixed    <NA>  TreatmentNuclei_25:Plant…   -0.199     0.226    -0.88    0.379
 8 fixed    <NA>  TreatmentNuclei_50:Plant…   -0.77      0.193    -3.99    0    
 9 fixed    <NA>  TreatmentStrip_25:Plante…   -1.46      0.205    -7.12    0    
10 ran_pars Block sd__(Intercept)              0.564    NA        NA      NA    
show code
sjPlot::tab_model(mod_binomial_full)
  cbind(n planted ind,n
ind-n planted ind)
Predictors Odds Ratios CI p
(Intercept) 0.10 0.06 – 0.17 <0.001
Treatment [Nuclei_25] 0.35 0.24 – 0.51 <0.001
Treatment [Nuclei_50] 0.75 0.56 – 1.00 0.051
Treatment [Strip_25] 0.90 0.67 – 1.22 0.506
Treatment [Strip_50] 0.72 0.53 – 0.97 0.033
Planted [Outside] 2.90 2.26 – 3.72 <0.001
Treatment [Nuclei_25] ×
Planted [Outside]
0.82 0.53 – 1.28 0.379
Treatment [Nuclei_50] ×
Planted [Outside]
0.46 0.32 – 0.68 <0.001
Treatment [Strip_25] ×
Planted [Outside]
0.23 0.16 – 0.35 <0.001
Random Effects
σ2 0.05
τ00 Block 0.32
ICC 0.86
N Block 5
Observations 178
Marginal R2 / Conditional R2 0.354 / 0.910

(ii) Not considering Baccharis + Vernonanthura

Figure species richness

Hill q = 0

Only small regenerants and considering ruderal species (Vernonanthura & Baccharis)

show code
dados_list <- species_table |>
  dplyr::filter(!scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes"))|>
  dplyr::group_by(Treatment, Planted, scientificName) |>
  dplyr::summarise(n_total = sum(n_sapling), .groups = "drop") |>
  dplyr::group_by(Treatment, Planted) |>
  dplyr::summarise(abundancias = list(n_total), .groups = "drop") |>
  dplyr::mutate(nome_grupo = paste(Treatment, Planted, sep = "_")) %>% 
  { setNames(.$abundancias, .$nome_grupo) }
res_inext <- iNEXT::iNEXT(
  dados_list,
  q = 0,
  size=seq(0, 450, 10),
  datatype = "abundance"
)

estimativas <- res_inext$iNextEst |>
  dplyr::bind_rows() |>
  tidyr::separate(Assemblage, into = c("Treatment", "Percent", "Planted"), sep = "_")

curvas_df0 <- estimativas |> 
  dplyr::filter(Order.q==0)|> 
  dplyr::filter(!is.na(SC.LCL))|>
  dplyr::mutate(
    Percent = ifelse(Percent=="Planted", "100", Percent),
    Percent = ordered(Percent, levels = c("25", "50", "100")),
    Planted = ifelse(Planted=="Inside", "Planted", "Unplanted"),
    Treatment = ifelse(Treatment=="Full", "Planting 100%", Treatment),
    Treatment = factor(Treatment, levels = c("Nuclei", "Strip", "Planting 100%"), ordered = TRUE)
  ) 

fig_recruit_rich <- curvas_df0 |>
  ggplot(aes(x = m, y = qD, color = Treatment, linetype = Planted)) +
  geom_line(linewidth = 1) +
  geom_ribbon(aes(ymin = qD.LCL, ymax = qD.UCL, fill = Treatment), alpha = 0.2, colour = NA) +
   geom_point(
    data = dplyr::filter(curvas_df0, Method == "Observed"),
    size = 2.5
  ) +
  facet_grid(Percent~., 
             labeller =  as_labeller(
               c("25" = "(a) Planting 25%",
                 "50" = "(b) Planting 50%",
                 "100" = "(c) Planting 100%")
               )) +
  scale_color_manual(values=c("red", "blue", "darkgreen"))+
  scale_fill_manual(values=c("red", "blue", "darkgreen"))+
  labs(
    x = "Number of sampling units (regenerants)",
    y = "Species diversity (Hill Index, q=0)"
    ) +
  ggthemes::theme_clean()+
  scale_x_continuous(expand = c(0.001, 0), breaks = seq(0, 450, 50)) +
  scale_y_continuous(expand = c(0, 0))+
  theme(
     panel.grid.major.y = element_blank(),
    panel.grid.major.x = element_blank(),
    axis.text.x = element_text(color = "black", hjust = 1, size = 10),
    axis.text.y = element_text(color = "black", size = 9),
    axis.title.y = element_text(color = "black", size = 11),
    strip.text = element_text(color = "black", hjust = 0.5, size = 10),
    strip.placement = "outside",
    legend.text = element_text(size = 11),
    legend.title = element_blank(),
    legend.key.size = unit(5, "mm"),
    legend.box.background = element_blank(),
    legend.box.margin = margin(t = 5, r = 0, b = 0, l = 0, "mm"),
    legend.background = element_blank(),
    legend.margin = margin(0, 0, 0, 0),
    legend.position = c(0.25, 0.25),
    panel.border = element_rect(colour = "black", fill = NA),
    panel.spacing = unit(0, "lines")
  )

Plotting and saving

show code
fig_recruit_rich

show code
ggsave("fig_richness_no_ruderal_small_recruits.png", plot=fig_recruit_rich, width=90, height = 210 ,unit="mm")

Table S2 - Recruitment (no ruderal)

Estimating Density of recruits (\(ind.ha^{-1}\)) without Baccharis and Vernonanthura

show code
density_recruits_n_m2 <- species_table |> 
  dplyr::filter(!scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes"))|>
  dplyr::mutate(
    n_small = dplyr::if_else(is.na(n_sapling), 0, n_sapling),
    n_large = dplyr::if_else(is.na(n_tree), 0, n_tree)) |>
  dplyr::group_by(Treatment, Block, Planted) |>
  dplyr::summarise(
    n_small_m2 = (sum(n_small, na.rm=T))/(144*4), 
    n_large_m2 = (sum(n_large, na.rm=T))/(144*4), 
    .groups = "drop")

dens_ha <- density_recruits_n_m2 |>
  dplyr::mutate(
    n_small_ha = n_small_m2 * 10000,
    n_large_ha = n_large_m2 * 10000
  )
dens_P_UP <- dens_ha |>
  tidyr::pivot_wider(
    names_from = Planted,
    values_from = c(n_small_ha, n_large_ha),
    names_sep = "_"
  )

prop_table <- tidyr::tibble(
  Treatment = unique(sort(species_table$Treatment)),
  prop_P = c(1, 0.25, 0.5, 0.25, 0.5),
  prop_UP = 1 - prop_P
)

result_regenerant <- dens_P_UP |>
  dplyr::left_join(prop_table, by = "Treatment") |>
  dplyr::mutate(
    n_small_total_ha = dplyr::coalesce(n_small_ha_Inside * prop_P, 0)  + dplyr::coalesce(n_small_ha_Outside * prop_UP, 0),
    n_large_total_ha = dplyr::coalesce(n_large_ha_Inside * prop_P, 0) + dplyr::coalesce(n_large_ha_Outside * prop_UP, 0)
  ) |>
  dplyr::select(Treatment,
         starts_with("n_small_ha"),
         starts_with("n_large_ha"),
         n_small_total_ha,
         n_large_total_ha
  )

Smaller recruits

  • Only in planted areas (without ruderal species)
show code
result_regenerant |>
  dplyr::group_by(Treatment)|>
  dplyr::summarise(
    mean = round(mean(n_small_ha_Inside, na.rm=TRUE),1),
    min = round(min(n_small_ha_Inside, na.rm=TRUE),0),
    max = round(max(n_small_ha_Inside, na.rm=TRUE),0)
  )|>
  dplyr::mutate(
    larger = sprintf("%.1f (%.1f - %.1f)", 
                                mean, 
                                min, 
                                max)) |>
  dplyr::select(Treatment, larger) |>
  tidyr::pivot_wider(
    names_from = Treatment,
    values_from = larger
    )
# A tibble: 1 × 5
  Full_Planted           Nuclei_25             Nuclei_50       Strip_25 Strip_50
  <chr>                  <chr>                 <chr>           <chr>    <chr>   
1 503.5 (139.0 - 1146.0) 482.6 (174.0 - 903.0) 961.8 (278.0 -… 722.2 (… 614.6 (…
  • Only in unplanted areas (without ruderal species)
show code
result_regenerant |>
  dplyr::group_by(Treatment)|>
  dplyr::summarise(
    mean = round(mean(n_small_ha_Outside, na.rm=TRUE), 1),
    min = round(min(n_small_ha_Outside, na.rm=TRUE), 0),
    max = round(max(n_small_ha_Outside, na.rm=TRUE), 0)
  )|>
  dplyr::mutate(
    larger = sprintf("%.1f (%.1f - %.1f)", 
                                mean, 
                                min, 
                                max)) |>
  dplyr::select(Treatment, larger) |>
  tidyr::pivot_wider(
    names_from = Treatment,
    values_from = larger
    )
Warning: There were 2 warnings in `dplyr::summarise()`.
The first warning was:
ℹ In argument: `min = round(min(n_small_ha_Outside, na.rm = TRUE), 0)`.
ℹ In group 1: `Treatment = "Full_Planted"`.
Caused by warning in `min()`:
! nenhum argumento não faltante para min; retornando Inf
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
# A tibble: 1 × 5
  Full_Planted     Nuclei_25               Nuclei_50           Strip_25 Strip_50
  <chr>            <chr>                   <chr>               <chr>    <chr>   
1 NaN (Inf - -Inf) 1177.1 (469.0 - 1910.0) 822.9 (538.0 - 138… 1145.8 … 1381.9 …
  • In total (planted + unplanted areas)
show code
result_regenerant |>
  dplyr::group_by(Treatment)|>
  dplyr::summarise(
    mean = round(mean(n_small_total_ha, na.rm=TRUE),1),
    min = round(min(n_small_total_ha, na.rm=TRUE),0),
    max = round(max(n_small_total_ha, na.rm=TRUE),0)
  )|>
  dplyr::mutate(
    larger = sprintf("%.1f (%.1f - %.1f)", 
                                mean, 
                                min, 
                                max)) |>
  dplyr::select(Treatment, larger) |>
  tidyr::pivot_wider(
    names_from = Treatment,
    values_from = larger
    )
# A tibble: 1 × 5
  Full_Planted           Nuclei_25             Nuclei_50       Strip_25 Strip_50
  <chr>                  <chr>                 <chr>           <chr>    <chr>   
1 503.5 (139.0 - 1146.0) 501.7 (43.0 - 1432.0) 446.2 (139.0 -… 577.7 (… 499.1 (…

Larger recruits

  • Only in planted areas (without ruderal species)
show code
result_regenerant |>
  dplyr::group_by(Treatment)|>
  dplyr::summarise(
    mean = round(mean(n_large_ha_Inside, na.rm=TRUE),1),
    min = round(min(n_large_ha_Inside, na.rm=TRUE),0),
    max = round(max(n_large_ha_Inside, na.rm=TRUE),0)
  ) |>
  dplyr::mutate(
    larger = sprintf("%.1f (%.1f - %.1f)", 
                                mean, 
                                min, 
                                max)) |>
  dplyr::select(Treatment, larger) |>
  tidyr::pivot_wider(
    names_from = Treatment,
    values_from = larger
    )
# A tibble: 1 × 5
  Full_Planted         Nuclei_25           Nuclei_50           Strip_25 Strip_50
  <chr>                <chr>               <chr>               <chr>    <chr>   
1 208.3 (87.0 - 486.0) 211.8 (0.0 - 608.0) 159.7 (35.0 - 278.… 149.3 (… 222.2 (…
  • Only in unplanted areas
show code
result_regenerant |>
  dplyr::group_by(Treatment)|>
  dplyr::summarise(
    mean = round(mean(n_large_ha_Outside, na.rm=TRUE),1),
    min = round(min(n_large_ha_Outside, na.rm=TRUE),0),
    max = round(max(n_large_ha_Outside, na.rm=TRUE),0)
  )|>
  dplyr::mutate(
    larger = sprintf("%.1f (%.1f - %.1f)", 
                                mean, 
                                min, 
                                max)) |>
  dplyr::select(Treatment, larger) |>
  tidyr::pivot_wider(
    names_from = Treatment,
    values_from = larger
    )
Warning: There were 2 warnings in `dplyr::summarise()`.
The first warning was:
ℹ In argument: `min = round(min(n_large_ha_Outside, na.rm = TRUE), 0)`.
ℹ In group 1: `Treatment = "Full_Planted"`.
Caused by warning in `min()`:
! nenhum argumento não faltante para min; retornando Inf
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
# A tibble: 1 × 5
  Full_Planted     Nuclei_25           Nuclei_50             Strip_25   Strip_50
  <chr>            <chr>               <chr>                 <chr>      <chr>   
1 NaN (Inf - -Inf) 333.3 (0.0 - 851.0) 138.9 (104.0 - 208.0) 253.5 (35… 364.6 (…
  • In total (planted + unplanted areas)
show code
result_regenerant |>
  dplyr::group_by(Treatment)|>
  dplyr::summarise(
    mean = round(mean(n_large_total_ha, na.rm = TRUE), 1),
    min = round(min(n_large_total_ha, na.rm = TRUE), 0),
    max = round(max(n_large_total_ha, na.rm = TRUE), 0)
  )|>
  dplyr::mutate(
    larger = sprintf("%.1f (%.1f - %.1f)", 
                                mean, 
                                min, 
                                max)) |>
  dplyr::select(Treatment, larger) |>
  tidyr::pivot_wider(
    names_from = Treatment,
    values_from = larger
    )
# A tibble: 1 × 5
  Full_Planted         Nuclei_25           Nuclei_50           Strip_25 Strip_50
  <chr>                <chr>               <chr>               <chr>    <chr>   
1 208.3 (87.0 - 486.0) 151.5 (0.0 - 638.0) 74.7 (17.0 - 139.0) 126.4 (… 146.7 (…

Analysis for smaller recruits abundance

Data:

show code
small_df <- species_table |>
  dplyr::filter(!scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes"))|>
  dplyr::group_by(Treatment, Block, Plot, Planted)|>
  dplyr::summarise(n = sum(n_sapling), .groups="drop")

small_df |>
  ggplot(aes(x=Planted, y=n, color=Block))+
  geom_point()+facet_wrap(~Treatment)+theme_bw()

Models:

show code
modelo_reg1 <- glmmTMB::glmmTMB(n ~ Treatment*Planted + (1 | Block),
                   # family = poisson(),
                   family = glmmTMB::nbinom2(),
                   data = small_df)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
show code
modelo_reg2 <-  glmmTMB::glmmTMB(n ~ Treatment + (1 | Block),
                   # family = poisson(),
                   family = glmmTMB::nbinom2(),
                   data = small_df)
modelo_reg3 <-  glmmTMB::glmmTMB(n ~ Planted + (1 | Block),
                   # family = poisson(),
                   family = glmmTMB::nbinom2(),
                   data = small_df)
modelo_reg_null <-  glmmTMB::glmmTMB(n ~ 0 + (1 | Block),
                   # family = poisson(), # Poisson
                   family = glmmTMB::nbinom2(), # NB
                   data = small_df)

# NBinom 1270.77 (null) |   1226.20 (Treatment*Planted) | 1250.26   (only Treament) | 1227.00   (only Planted)
# Poisson  4032.177 |   3373.854        |   3813.769        |   3759.763

AIC(modelo_reg_null, modelo_reg1, modelo_reg2, modelo_reg3)
                df      AIC
modelo_reg_null  2 1270.773
modelo_reg1     11 1226.199
modelo_reg2      7 1250.263
modelo_reg3      4 1227.001

non-positive-definite Hessian matrix indicates that there are problems in estimating the model parameters. This may be related to redundancies in the random effects, perfect separation, few data for some combinations of factors, or problems in choosing the distribution. Therefore, let’s check the singularity:

show code
performance::check_singularity(modelo_reg1)
[1] FALSE
show code
performance::check_model(modelo_reg1)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
`check_outliers()` does not yet support models of class `glmmTMB`.

Exploring residuals

show code
par(mfrow=c(2,2))
DHARMa::plotResiduals(modelo_reg1)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
show code
DHARMa::plotQQunif(modelo_reg1)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
show code
DHARMa::plotResiduals(modelo_reg3)
DHARMa::plotQQunif(modelo_reg3)

show code
car::Anova(modelo_reg1)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: n
                    Chisq Df Pr(>Chisq)    
Treatment          4.7099  4    0.31838    
Planted           24.8150  1  6.311e-07 ***
Treatment:Planted 10.7731  3    0.01302 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Is there any difference between treatments considering only inside or outside?

show code
emmeans::emmeans(modelo_reg1, ~ Treatment|Planted) |> 
  multcomp::cld(adjust="bonferroni")
Planted = Inside:
 Treatment    emmean    SE  df asymp.LCL asymp.UCL .group
 Full_Planted   1.86 0.241 Inf      1.24      2.48  1    
 Nuclei_25      1.87 0.241 Inf      1.25      2.50  1    
 Strip_50       2.09 0.238 Inf      1.48      2.70  1    
 Strip_25       2.27 0.236 Inf      1.66      2.88  1    
 Nuclei_50      2.42 0.234 Inf      1.82      3.02  1    

Planted = Outside:
 Treatment    emmean    SE  df asymp.LCL asymp.UCL .group
 Nuclei_50      2.43 0.234 Inf      1.84      3.01  1    
 Strip_25       2.73 0.232 Inf      2.15      3.31  1    
 Nuclei_25      2.79 0.231 Inf      2.21      3.36  1    
 Strip_50       2.99 0.230 Inf      2.41      3.56  1    
 Full_Planted nonEst    NA  NA        NA        NA       

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for varying numbers of estimates 
P value adjustment: bonferroni method for varying numbers of tests 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Answer: No

Is there any difference between inside/outside considering treatments?

show code
emmeans::emmeans(modelo_reg1, ~ Planted|Treatment) |> 
    multcomp::cld(adjust="bonferroni")
Treatment = Full_Planted:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Inside    1.86 0.241 Inf      1.39      2.33  1    
 Outside nonEst    NA  NA        NA        NA       

Treatment = Nuclei_25:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Inside    1.87 0.241 Inf      1.33      2.42  1    
 Outside   2.79 0.231 Inf      2.27      3.30   2   

Treatment = Nuclei_50:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Inside    2.42 0.234 Inf      1.90      2.95  1    
 Outside   2.43 0.234 Inf      1.90      2.95  1    

Treatment = Strip_25:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Inside    2.27 0.236 Inf      1.74      2.80  1    
 Outside   2.73 0.232 Inf      2.21      3.25   2   

Treatment = Strip_50:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Inside    2.09 0.238 Inf      1.56      2.63  1    
 Outside   2.99 0.230 Inf      2.47      3.51   2   

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for varying numbers of estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Answer: Yes.

show code
sjPlot::tab_model(modelo_reg1)
Model matrix is rank deficient. Parameters
  `TreatmentStrip_50:PlantedOutside` were not estimable.
  n
Predictors Incidence Rate Ratios CI p
(Intercept) 6.43 4.01 – 10.31 <0.001
Treatment [Nuclei_25] 1.01 0.63 – 1.62 0.955
Treatment [Nuclei_50] 1.75 1.11 – 2.76 0.016
Treatment [Strip_25] 1.51 0.95 – 2.39 0.079
Treatment [Strip_50] 1.26 0.79 – 2.00 0.326
Planted [Outside] 2.45 1.57 – 3.82 <0.001
Treatment [Nuclei_25] ×
Planted [Outside]
1.01 0.54 – 1.91 0.964
Treatment [Nuclei_50] ×
Planted [Outside]
0.41 0.22 – 0.77 0.005
Treatment [Strip_25] ×
Planted [Outside]
0.65 0.35 – 1.21 0.171
Random Effects
σ2 0.40
τ00 Block 0.15
ICC 0.27
N Block 5
Observations 180
Marginal R2 / Conditional R2 0.207 / 0.419
show code
broom.mixed::tidy(modelo_reg1) |>
  dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3))) |> dplyr::select(-c(group,component))
# A tibble: 11 × 6
   effect   term                            estimate std.error statistic p.value
   <chr>    <chr>                              <dbl>     <dbl>     <dbl>   <dbl>
 1 fixed    (Intercept)                        1.86      0.241     7.72    0    
 2 fixed    TreatmentNuclei_25                 0.014     0.24      0.056   0.955
 3 fixed    TreatmentNuclei_50                 0.56      0.232     2.41    0.016
 4 fixed    TreatmentStrip_25                  0.411     0.234     1.75    0.079
 5 fixed    TreatmentStrip_50                  0.232     0.236     0.982   0.326
 6 fixed    PlantedOutside                     0.897     0.226     3.97    0    
 7 fixed    TreatmentNuclei_25:PlantedOuts…    0.015     0.322     0.045   0.964
 8 fixed    TreatmentNuclei_50:PlantedOuts…   -0.892     0.319    -2.80    0.005
 9 fixed    TreatmentStrip_25:PlantedOutsi…   -0.437     0.319    -1.37    0.171
10 fixed    TreatmentStrip_50:PlantedOutsi…   NA        NA        NA      NA    
11 ran_pars sd__(Intercept)                    0.383    NA        NA      NA    

Analysis for larger recruits abundance

Data:

show code
large_regenerants_2024 <- species_table |>
  dplyr::filter(!scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes"))|>
  dplyr::group_by(Treatment, Block, Plot, Planted)|>
  dplyr::summarise(n = sum(n_tree), .groups="drop")

Plotting

show code
large_regenerants_2024 |>
  ggplot(aes(x=Planted, y=n, color=Block))+
  geom_point()+facet_wrap(~Treatment)+theme_bw()

Models

show code
modelo_reg1 <- glmmTMB::glmmTMB(n ~ Treatment*Planted + (1 | Block),
                   ziformula = ~1,
                    # family = poisson(),
                     family = glmmTMB::nbinom2(),
                    # family = glmmTMB::truncated_nbinom2(),
                   data = large_regenerants_2024)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
show code
modelo_reg2 <-  glmmTMB::glmmTMB(n ~ Treatment + (1 | Block),
                   ziformula = ~1,
                   #  family = poisson(),
                    family = glmmTMB::nbinom2(),
                   # family = glmmTMB::truncated_nbinom2(),
                   data = large_regenerants_2024)
modelo_reg3 <-  glmmTMB::glmmTMB(n ~ Planted + (1 | Block),
                   ziformula = ~1,
                   # family = poisson(),
                    family = glmmTMB::nbinom2(),
                   #  family = glmmTMB::truncated_nbinom2(),
                   data = large_regenerants_2024)
modelo_reg_null <-  glmmTMB::glmmTMB(n ~ 0 + (1 | Block),
                   ziformula = ~1,
                    # family = poisson(), # ZIP
                     family = glmmTMB::nbinom2(), # ZINB
                    # family = glmmTMB::truncated_nbinom2(), # hurdle
                   data = large_regenerants_2024)
# Model (null) | (treat*planted) | (treat) | (planted)
# ZINB 855.21       |   843.24       |  840.70  |   838.42

AIC(modelo_reg_null, modelo_reg1, modelo_reg2, modelo_reg3)
                df      AIC
modelo_reg_null  3 855.2141
modelo_reg1     12 843.2460
modelo_reg2      8 840.7029
modelo_reg3      5 838.4200
show code
performance::check_model(modelo_reg1)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
`check_outliers()` does not yet support models of class `glmmTMB`.

Exploring residuals

show code
par(mfrow=c(2,2))
DHARMa::plotResiduals(modelo_reg1)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
show code
DHARMa::plotQQunif(modelo_reg1)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
show code
DHARMa::plotResiduals(modelo_reg3)
DHARMa::plotQQunif(modelo_reg3)

show code
car::Anova(modelo_reg1)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: n
                   Chisq Df Pr(>Chisq)  
Treatment         7.9554  4    0.09323 .
Planted           3.4787  1    0.06216 .
Treatment:Planted 2.1403  3    0.54381  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Is there any difference between treatments considering inside/outside?

show code
emmeans::emmeans(modelo_reg1, ~ Treatment|Planted) |> 
  multcomp::cld(adjust = "bonferroni")
Planted = Inside:
 Treatment    emmean    SE  df asymp.LCL asymp.UCL .group
 Strip_25      0.920 0.292 Inf     0.170      1.67  1    
 Nuclei_50     0.971 0.282 Inf     0.245      1.70  1    
 Full_Planted  1.162 0.249 Inf     0.520      1.80  1    
 Strip_50      1.278 0.261 Inf     0.606      1.95  1    
 Nuclei_25     1.280 0.280 Inf     0.558      2.00  1    

Planted = Outside:
 Treatment    emmean    SE  df asymp.LCL asymp.UCL .group
 Nuclei_50     0.820 0.284 Inf     0.109      1.53  1    
 Strip_25      1.419 0.259 Inf     0.772      2.07  1    
 Strip_50      1.744 0.238 Inf     1.150      2.34  1    
 Nuclei_25     1.787 0.283 Inf     1.079      2.49  1    
 Full_Planted nonEst    NA  NA        NA        NA       

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for varying numbers of estimates 
P value adjustment: bonferroni method for varying numbers of tests 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Answer: No.

Is there any difference between inside/outside considering treatments?

show code
emmeans::emmeans(modelo_reg1, ~ Planted|Treatment) |> 
  multcomp::cld(adjust = "bonferroni")
Treatment = Full_Planted:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Inside   1.162 0.249 Inf     0.674      1.65  1    
 Outside nonEst    NA  NA        NA        NA       

Treatment = Nuclei_25:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Inside   1.280 0.280 Inf     0.652      1.91  1    
 Outside  1.787 0.283 Inf     1.152      2.42  1    

Treatment = Nuclei_50:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Outside  0.820 0.284 Inf     0.182      1.46  1    
 Inside   0.971 0.282 Inf     0.339      1.60  1    

Treatment = Strip_25:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Inside   0.920 0.292 Inf     0.267      1.57  1    
 Outside  1.419 0.259 Inf     0.838      2.00  1    

Treatment = Strip_50:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Inside   1.278 0.261 Inf     0.694      1.86  1    
 Outside  1.744 0.238 Inf     1.211      2.28  1    

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for varying numbers of estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Answer: No.

Table S3

show code
sjPlot::tab_model(modelo_reg1)
Model matrix is rank deficient. Parameters
  `TreatmentStrip_50:PlantedOutside` were not estimable.
  n
Predictors Incidence Rate Ratios CI p
Count Model
(Intercept) 3.20 1.96 – 5.21 <0.001
Treatment [Nuclei_25] 1.13 0.55 – 2.31 0.747
Treatment [Nuclei_50] 0.83 0.40 – 1.70 0.605
Treatment [Strip_25] 0.79 0.38 – 1.64 0.520
Treatment [Strip_50] 1.12 0.56 – 2.25 0.743
Planted [Outside] 1.59 0.81 – 3.14 0.178
Treatment [Nuclei_25] ×
Planted [Outside]
1.04 0.38 – 2.82 0.937
Treatment [Nuclei_50] ×
Planted [Outside]
0.54 0.20 – 1.49 0.233
Treatment [Strip_25] ×
Planted [Outside]
1.03 0.38 – 2.80 0.949
(Intercept) 3.53 2.03 – 9.52
Zero-Inflated Model
(Intercept) 0.14 0.04 – 0.50 0.003
Random Effects
σ2
τ00 Block 0.00
N Block 5
Observations 180
show code
broom.mixed::tidy(modelo_reg1) |>
  dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3))) |>
  knitr::kable()
effect component group term estimate std.error statistic p.value
fixed cond NA (Intercept) 1.162 0.249 4.663 0.000
fixed cond NA TreatmentNuclei_25 0.119 0.367 0.323 0.747
fixed cond NA TreatmentNuclei_50 -0.191 0.369 -0.517 0.605
fixed cond NA TreatmentStrip_25 -0.241 0.376 -0.643 0.520
fixed cond NA TreatmentStrip_50 0.116 0.355 0.328 0.743
fixed cond NA PlantedOutside 0.466 0.346 1.347 0.178
fixed cond NA TreatmentNuclei_25:PlantedOutside 0.040 0.508 0.080 0.937
fixed cond NA TreatmentNuclei_50:PlantedOutside -0.617 0.518 -1.192 0.233
fixed cond NA TreatmentStrip_25:PlantedOutside 0.033 0.508 0.064 0.949
fixed cond NA TreatmentStrip_50:PlantedOutside NA NA NA NA
fixed zi NA (Intercept) -1.962 0.652 -3.007 0.003
ran_pars cond Block sd__(Intercept) 0.000 NA NA NA

Propotion of smaller recruits from planted species

  • Without ruderal species

Getting proportions of planted/unplanted species

show code
prop_regenerant_part <- species_table |>
  dplyr::filter(!scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes"))|>
  dplyr::mutate(
    planted_tree = ifelse(`Planting group`=="NA", "No", "Yes"))|>
  dplyr::filter(n_sapling>0)|>
  dplyr::group_by(Treatment, Planted, Block, Plot) |>
  dplyr::summarise(
    n_species = dplyr::n(),
    n_ind = sum(n_sapling),
    n_planted_ind = sum(ifelse(planted_tree == "Yes", n_sapling, 0), na.rm = TRUE),
    n_planted_ind = ifelse(is.na(n_planted_ind), 0, n_planted_ind),
    n_planted_spp = sum(planted_tree == "Yes"), 
    n_planted_spp = ifelse(is.na(n_planted_spp), 0, n_planted_spp),
    .groups = "drop") |>
  dplyr::mutate(
    prop_planted_ind = n_planted_ind / n_ind,
    prop_planted_spp = n_planted_spp / n_species
  )

prop_regenerant_part |>
  dplyr::group_by(Treatment, Planted)|>
  dplyr::summarise(
    prop_planted_mean = mean(prop_planted_ind, na.rm=T),
    prop_planted_min = min(prop_planted_ind, na.rm=T),
    prop_planted_max = max(prop_planted_ind, na.rm=T),
    prop_planted_spp_mean = mean(prop_planted_spp, na.rm=T),
    prop_planted_spp_min = min(prop_planted_spp, na.rm=T),
    prop_planted_spp_max = max(prop_planted_spp, na.rm=T),
    .groups = "drop"
  ) |>
  dplyr::mutate(
    prop_planted = sprintf("%.1f (%.1f - %.1f)", 
                                round(prop_planted_mean,3)*100, 
                                round(prop_planted_min,3)*100, 
                                round(prop_planted_max,3)*100),
    prop_planted_spp = sprintf("%.1f (%.1f - %.1f)", 
                                    round(prop_planted_spp_mean,3)*100, 
                                    round(prop_planted_spp_min,3)*100, 
                                    round(prop_planted_spp_max,3)*100)
    ) |>
  dplyr::select(Treatment,Planted, prop_planted, prop_planted_spp) |>
  tidyr::pivot_wider(
    names_from = Treatment,
    values_from = c(prop_planted, prop_planted_spp)
    ) # |> writexl::write_xlsx("resulta_prop_recruit.xlsx")
# A tibble: 2 × 11
  Planted prop_planted_Full_Plan…¹ prop_planted_Nuclei_25 prop_planted_Nuclei_50
  <chr>   <chr>                    <chr>                  <chr>                 
1 Inside  54.4 (0.0 - 100.0)       41.4 (0.0 - 100.0)     37.5 (0.0 - 100.0)    
2 Outside <NA>                     25.0 (0.0 - 91.2)      27.0 (0.0 - 92.3)     
# ℹ abbreviated name: ¹​prop_planted_Full_Planted
# ℹ 7 more variables: prop_planted_Strip_25 <chr>, prop_planted_Strip_50 <chr>,
#   prop_planted_spp_Full_Planted <chr>, prop_planted_spp_Nuclei_25 <chr>,
#   prop_planted_spp_Nuclei_50 <chr>, prop_planted_spp_Strip_25 <chr>,
#   prop_planted_spp_Strip_50 <chr>
show code
prop_regenerant_part |>
  dplyr::group_by(Planted)|>
  dplyr::summarise(
    prop_planted_mean = mean(prop_planted_ind, na.rm=T),
    prop_planted_se = (sd(prop_planted_ind, na.rm=T))/sqrt(dplyr::n()),
    prop_planted_spp_mean = mean(prop_planted_spp, na.rm=T),
    prop_planted_spp_se = (sd(prop_planted_spp, na.rm=T))/sqrt(dplyr::n()),
    .groups = "drop"
  )
# A tibble: 2 × 5
  Planted prop_planted_mean prop_planted_se prop_planted_spp_mean
  <chr>               <dbl>           <dbl>                 <dbl>
1 Inside              0.524          0.0372                 0.488
2 Outside             0.283          0.0309                 0.243
# ℹ 1 more variable: prop_planted_spp_se <dbl>

Statistical analysis for proportion of individuals:

show code
mod_binomial_part <-lme4::glmer(
  cbind(n_planted_ind, n_ind - n_planted_ind) ~ Treatment*Planted + (1 | Block),
  family = binomial(link = "logit"),
  data = prop_regenerant_part)
fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
show code
performance::check_overdispersion(mod_binomial_part) # yes
# Overdispersion test

 dispersion ratio =   1.837
          p-value = < 0.001
Overdispersion detected.

It should be Beta

show code
mod_beta_part <-glmmTMB::glmmTMB(
  cbind(n_planted_ind, n_ind - n_planted_ind) ~ Treatment*Planted + (1 | Block),
  family = glmmTMB::betabinomial(link = "logit"),
  data = prop_regenerant_part)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
show code
performance::check_model(mod_beta_part)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
`check_outliers()` does not yet support models of class `glmmTMB`.

show code
DHARMa::plotResiduals(mod_beta_part)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside

Pairwise:

show code
emmeans::emmeans(mod_beta_part, pairwise ~ Treatment | Planted) |> multcomp::cld(adjust="sidak") # no differences
Planted = Inside:
 Treatment    emmean    SE  df asymp.LCL asymp.UCL .group
 Nuclei_50    -0.429 0.311 Inf    -1.228     0.370  1    
 Nuclei_25    -0.358 0.314 Inf    -1.165     0.449  1    
 Full_Planted  0.243 0.326 Inf    -0.595     1.082  1    
 Strip_25      0.447 0.298 Inf    -0.319     1.212  1    
 Strip_50      0.490 0.292 Inf    -0.260     1.240  1    

Planted = Outside:
 Treatment    emmean    SE  df asymp.LCL asymp.UCL .group
 Strip_25     -1.196 0.298 Inf    -1.938    -0.455  1    
 Nuclei_25    -1.096 0.287 Inf    -1.810    -0.382  1    
 Nuclei_50    -0.935 0.314 Inf    -1.718    -0.153  1    
 Strip_50     -0.376 0.264 Inf    -1.035     0.282  1    
 Full_Planted nonEst    NA  NA        NA        NA       

Results are given on the logit (not the response) scale. 
Confidence level used: 0.95 
Conf-level adjustment: sidak method for varying numbers of estimates 
Results are given on the log odds ratio (not the response) scale. 
P value adjustment: sidak method for varying numbers of tests 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
show code
emmeans::emmeans(mod_beta_part, pairwise ~  Planted | Treatment) |> multcomp::cld(adjust="sidak") # All different
Treatment = Full_Planted:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Inside   0.243 0.326 Inf    -0.397     0.883  1    
 Outside nonEst    NA  NA        NA        NA       

Treatment = Nuclei_25:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Outside -1.096 0.287 Inf    -1.737    -0.455  1    
 Inside  -0.358 0.314 Inf    -1.061     0.345  1    

Treatment = Nuclei_50:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Outside -0.935 0.314 Inf    -1.638    -0.233  1    
 Inside  -0.429 0.311 Inf    -1.125     0.267  1    

Treatment = Strip_25:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Outside -1.196 0.298 Inf    -1.862    -0.530  1    
 Inside   0.447 0.298 Inf    -0.220     1.113   2   

Treatment = Strip_50:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Outside -0.376 0.264 Inf    -0.968     0.215  1    
 Inside   0.490 0.292 Inf    -0.163     1.143   2   

Results are given on the logit (not the response) scale. 
Confidence level used: 0.95 
Conf-level adjustment: sidak method for varying numbers of estimates 
Results are given on the log odds ratio (not the response) scale. 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Statistical analysis for proportion of individuals:

show code
mod_binomial_prop <-lme4::glmer(
  cbind(n_planted_spp, n_species - n_planted_spp) ~ Treatment*Planted + (1 | Block),
  family = binomial(link = "logit"),
  data = prop_regenerant_part)
fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
boundary (singular) fit: see help('isSingular')
show code
performance::check_overdispersion(mod_binomial_prop) # underdispersion
# Overdispersion test

 dispersion ratio = 0.769
          p-value = 0.032
Underdispersion detected.

It should be Beta

show code
mod_beta_prop <-glmmTMB::glmmTMB(
  cbind(n_planted_spp, n_species - n_planted_spp) ~ Treatment*Planted + (1 | Block),
  family = glmmTMB::betabinomial(link = "logit"),
  data = prop_regenerant_part)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
problem; false convergence (8). See vignette('troubleshooting'),
help('diagnose')
show code
performance::check_model(mod_beta_prop)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
`check_outliers()` does not yet support models of class `glmmTMB`.

show code
DHARMa::plotResiduals(mod_beta_prop)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside

Pairwise:

show code
emmeans::emmeans(mod_beta_prop, pairwise ~ Treatment | Planted) |> multcomp::cld(adjust="sidak") # no differences
Planted = Inside:
 Treatment    emmean    SE  df asymp.LCL asymp.UCL .group
 Nuclei_25    -0.582 0.286 Inf    -1.318     0.154  1    
 Nuclei_50    -0.528 0.263 Inf    -1.203     0.147  1    
 Strip_25     -0.210 0.246 Inf    -0.841     0.421  1    
 Full_Planted  0.251 0.291 Inf    -0.496     0.999  1    
 Strip_50      0.297 0.259 Inf    -0.368     0.962  1    

Planted = Outside:
 Treatment    emmean    SE  df asymp.LCL asymp.UCL .group
 Nuclei_25    -1.399 0.250 Inf    -2.021    -0.777  1    
 Strip_25     -1.204 0.233 Inf    -1.784    -0.624  1    
 Nuclei_50    -0.944 0.257 Inf    -1.585    -0.304  1    
 Strip_50     -0.847 0.208 Inf    -1.366    -0.329  1    
 Full_Planted nonEst    NA  NA        NA        NA       

Results are given on the logit (not the response) scale. 
Confidence level used: 0.95 
Conf-level adjustment: sidak method for varying numbers of estimates 
Results are given on the log odds ratio (not the response) scale. 
P value adjustment: sidak method for varying numbers of tests 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
show code
emmeans::emmeans(mod_beta_prop, pairwise ~  Planted | Treatment) |> multcomp::cld(adjust="sidak") # All different
Treatment = Full_Planted:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Inside   0.251 0.291 Inf    -0.319    0.8216  1    
 Outside nonEst    NA  NA        NA        NA       

Treatment = Nuclei_25:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Outside -1.399 0.250 Inf    -1.957   -0.8403  1    
 Inside  -0.582 0.286 Inf    -1.222    0.0587   2   

Treatment = Nuclei_50:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Outside -0.944 0.257 Inf    -1.520   -0.3693  1    
 Inside  -0.528 0.263 Inf    -1.116    0.0599  1    

Treatment = Strip_25:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Outside -1.204 0.233 Inf    -1.724   -0.6835  1    
 Inside  -0.210 0.246 Inf    -0.759    0.3398   2   

Treatment = Strip_50:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Outside -0.847 0.208 Inf    -1.313   -0.3820  1    
 Inside   0.297 0.259 Inf    -0.282    0.8763   2   

Results are given on the logit (not the response) scale. 
Confidence level used: 0.95 
Conf-level adjustment: sidak method for varying numbers of estimates 
Results are given on the log odds ratio (not the response) scale. 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Table S3

show code
broom.mixed::tidy(mod_beta_part) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))
# A tibble: 11 × 8
   effect   component group term            estimate std.error statistic p.value
   <chr>    <chr>     <chr> <chr>              <dbl>     <dbl>     <dbl>   <dbl>
 1 fixed    cond      <NA>  (Intercept)        0.243     0.326     0.745   0.456
 2 fixed    cond      <NA>  TreatmentNucle…   -0.601     0.454    -1.32    0.185
 3 fixed    cond      <NA>  TreatmentNucle…   -0.672     0.451    -1.49    0.136
 4 fixed    cond      <NA>  TreatmentStrip…    0.203     0.442     0.461   0.645
 5 fixed    cond      <NA>  TreatmentStrip…    0.247     0.437     0.564   0.572
 6 fixed    cond      <NA>  PlantedOutside    -0.867     0.394    -2.20    0.028
 7 fixed    cond      <NA>  TreatmentNucle…    0.129     0.577     0.223   0.824
 8 fixed    cond      <NA>  TreatmentNucle…    0.36      0.592     0.608   0.543
 9 fixed    cond      <NA>  TreatmentStrip…   -0.776     0.576    -1.35    0.178
10 fixed    cond      <NA>  TreatmentStrip…   NA        NA        NA      NA    
11 ran_pars cond      Block sd__(Intercept)    0        NA        NA      NA    
show code
sjPlot::tab_model(mod_beta_part)
Model matrix is rank deficient. Parameters
  `TreatmentStrip_50:PlantedOutside` were not estimable.
  cbind(n planted ind,n
ind-n planted ind)
Predictors Odds Ratios CI p
(Intercept) 1.28 0.67 – 2.42 0.456
Treatment [Nuclei_25] 0.55 0.23 – 1.33 0.185
Treatment [Nuclei_50] 0.51 0.21 – 1.24 0.136
Treatment [Strip_25] 1.23 0.52 – 2.91 0.645
Treatment [Strip_50] 1.28 0.54 – 3.02 0.572
Planted [Outside] 0.42 0.19 – 0.91 0.028
Treatment [Nuclei_25] ×
Planted [Outside]
1.14 0.37 – 3.52 0.824
Treatment [Nuclei_50] ×
Planted [Outside]
1.43 0.45 – 4.57 0.543
Treatment [Strip_25] ×
Planted [Outside]
0.46 0.15 – 1.42 0.178
Random Effects
σ2 0.37
τ00 Block 0.00
N Block 5
Observations 176
Marginal R2 / Conditional R2 0.500 / NA

(iii) Only considering Baccharis + Vernonanthura

Table S2 - Only ruderal

Estimating density of recruits (\(ind.ha^{-1}\)) and without ruderal (Baccharis and Vernonanthura)

show code
density_ruderal_n_m2 <- species_table |> 
  dplyr::filter(scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes"))|>
  # dplyr::select(Treatment, Block, Plot, Planted)|>unique()|>   with(table(Block, Plot, by=Treatment)) # Full_Planted P2B5 & P3B1 (Inside)
  dplyr::add_row(
    scientificName = c(NA, NA),
    `Seed dispersal` = c(NA, NA),
    `Planting group` = c(NA, NA),
    Treatment = c("Full_Planted", "Full_Planted"),
    Block = c("B5", "B1"),
    Plot = c("P2", "P3"),
    Planted = c("Inside","Inside"),
    n_sapling= c(0,0),
    n_tree = c(0,0)
    ) |>
  dplyr::mutate(
    n_small = dplyr::if_else(is.na(n_sapling), 0, n_sapling),
    n_large = dplyr::if_else(is.na(n_tree), 0, n_tree)) |>
  dplyr::group_by(Treatment, Block, Planted) |>
  dplyr::summarise(
    n_small_m2 = (sum(n_small, na.rm=T))/(144*4), 
    n_large_m2 = (sum(n_large, na.rm=T))/(144*4), 
    .groups = "drop")

dens_ruderal_ha <- density_ruderal_n_m2 |>
  dplyr::mutate(
    n_small_ha = n_small_m2 * 10000,
    n_large_ha = n_large_m2 * 10000
  )
dens_ruderal_P_UP <- dens_ruderal_ha |>
  tidyr::pivot_wider(
    names_from = Planted,
    values_from = c(n_small_ha, n_large_ha),
    names_sep = "_"
  )

prop_table <- tidyr::tibble(
  Treatment = unique(sort(species_table$Treatment)),
  prop_P = c(1, 0.25, 0.5, 0.25, 0.5),
  prop_UP = 1 - prop_P
)

result_ruderal <- dens_ruderal_P_UP |>
  dplyr::left_join(prop_table, by = "Treatment") |>
  dplyr::mutate(
    n_small_total_ha = dplyr::coalesce(n_small_ha_Inside * prop_P, 0)  + dplyr::coalesce(n_small_ha_Outside * prop_UP, 0),
    n_large_total_ha = dplyr::coalesce(n_large_ha_Inside * prop_P, 0) + dplyr::coalesce(n_large_ha_Outside * prop_UP, 0)
  ) |>
  dplyr::select(Treatment,
         starts_with("n_small_ha"),
         starts_with("n_large_ha"),
         n_small_total_ha,
         n_large_total_ha
  )

Smaller recruits

  • Only in planted areas
show code
result_ruderal |>
  dplyr::group_by(Treatment)|>
  dplyr::summarise(
    mean = round(mean(n_small_ha_Inside, na.rm=TRUE),1),
    min = round(min(n_small_ha_Inside, na.rm=TRUE),0),
    max = round(max(n_small_ha_Inside, na.rm=TRUE),0)
  )|>
  dplyr::mutate(
    larger = sprintf("%.1f (%.1f - %.1f)", 
                                mean, 
                                min, 
                                max)) |>
  dplyr::select(Treatment, larger) |>
  tidyr::pivot_wider(
    names_from = Treatment,
    values_from = larger
    )
# A tibble: 1 × 5
  Full_Planted             Nuclei_25                Nuclei_50  Strip_25 Strip_50
  <chr>                    <chr>                    <chr>      <chr>    <chr>   
1 2545.1 (1580.0 - 4566.0) 4388.9 (2847.0 - 6858.0) 5145.8 (1… 4093.8 … 5409.7 …
  • Only in unplanted areas
show code
result_ruderal |>
  dplyr::group_by(Treatment)|>
  dplyr::summarise(
    mean = round(mean(n_small_ha_Outside, na.rm=TRUE), 1),
    min = round(min(n_small_ha_Outside, na.rm=TRUE), 0),
    max = round(max(n_small_ha_Outside, na.rm=TRUE), 0)
  )|>
  dplyr::mutate(
    larger = sprintf("%.1f (%.1f - %.1f)", 
                                mean, 
                                min, 
                                max)) |>
  dplyr::select(Treatment, larger) |>
  tidyr::pivot_wider(
    names_from = Treatment,
    values_from = larger
    )
Warning: There were 2 warnings in `dplyr::summarise()`.
The first warning was:
ℹ In argument: `min = round(min(n_small_ha_Outside, na.rm = TRUE), 0)`.
ℹ In group 1: `Treatment = "Full_Planted"`.
Caused by warning in `min()`:
! nenhum argumento não faltante para min; retornando Inf
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
# A tibble: 1 × 5
  Full_Planted     Nuclei_25                Nuclei_50          Strip_25 Strip_50
  <chr>            <chr>                    <chr>              <chr>    <chr>   
1 NaN (Inf - -Inf) 2784.7 (1181.0 - 6858.0) 2638.9 (1302.0 - … 2246.5 … 2489.6 …
  • In total (planted + unplanted areas)
show code
result_ruderal |>
  dplyr::group_by(Treatment)|>
  dplyr::summarise(
    mean = round(mean(n_small_total_ha, na.rm=TRUE),1),
    min = round(min(n_small_total_ha, na.rm=TRUE),0),
    max = round(max(n_small_total_ha, na.rm=TRUE),0)
  )|>
  dplyr::mutate(
    larger = sprintf("%.1f (%.1f - %.1f)", 
                                mean, 
                                min, 
                                max)) |>
  dplyr::select(Treatment, larger) |>
  tidyr::pivot_wider(
    names_from = Treatment,
    values_from = larger
    )
# A tibble: 1 × 5
  Full_Planted             Nuclei_25               Nuclei_50   Strip_25 Strip_50
  <chr>                    <chr>                   <chr>       <chr>    <chr>   
1 2545.1 (1580.0 - 4566.0) 1592.9 (712.0 - 5143.0) 1946.2 (65… 1354.2 … 1974.8 …

Larger recruits

All treatments and inside/outside areas showed similar values

  • Only in planted areas
show code
result_ruderal |>
  dplyr::group_by(Treatment)|>
  dplyr::summarise(
    mean = round(mean(n_large_ha_Inside, na.rm=TRUE),1),
    min = round(min(n_large_ha_Inside, na.rm=TRUE),0),
    max = round(max(n_large_ha_Inside, na.rm=TRUE),0)
  ) |>
  dplyr::mutate(
    larger = sprintf("%.1f (%.1f - %.1f)", 
                                mean, 
                                min, 
                                max)) |>
  dplyr::select(Treatment, larger) |>
  tidyr::pivot_wider(
    names_from = Treatment,
    values_from = larger
    )
# A tibble: 1 × 5
  Full_Planted      Nuclei_25         Nuclei_50        Strip_25         Strip_50
  <chr>             <chr>             <chr>            <chr>            <chr>   
1 13.9 (0.0 - 35.0) 17.4 (0.0 - 69.0) 3.5 (0.0 - 17.0) 3.5 (0.0 - 17.0) 3.5 (0.…
  • Only in unplanted areas
show code
result_ruderal |>
  dplyr::group_by(Treatment)|>
  dplyr::summarise(
    mean = round(mean(n_large_ha_Outside, na.rm=TRUE),1),
    min = round(min(n_large_ha_Outside, na.rm=TRUE),0),
    max = round(max(n_large_ha_Outside, na.rm=TRUE),0)
  )|>
  dplyr::mutate(
    larger = sprintf("%.1f (%.1f - %.1f)", 
                                mean, 
                                min, 
                                max)) |>
  dplyr::select(Treatment, larger) |>
  tidyr::pivot_wider(
    names_from = Treatment,
    values_from = larger
    )
Warning: There were 2 warnings in `dplyr::summarise()`.
The first warning was:
ℹ In argument: `min = round(min(n_large_ha_Outside, na.rm = TRUE), 0)`.
ℹ In group 1: `Treatment = "Full_Planted"`.
Caused by warning in `min()`:
! nenhum argumento não faltante para min; retornando Inf
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
# A tibble: 1 × 5
  Full_Planted     Nuclei_25         Nuclei_50         Strip_25         Strip_50
  <chr>            <chr>             <chr>             <chr>            <chr>   
1 NaN (Inf - -Inf) 31.2 (0.0 - 52.0) 24.3 (0.0 - 87.0) 48.6 (0.0 - 139… 17.4 (0…
  • In total (planted + unplanted areas)
show code
result_ruderal |>
  dplyr::group_by(Treatment)|>
  dplyr::summarise(
    mean = round(mean(n_large_total_ha, na.rm = TRUE), 1),
    min = round(min(n_large_total_ha, na.rm = TRUE), 0),
    max = round(max(n_large_total_ha, na.rm = TRUE), 0)
  )|>
  dplyr::mutate(
    larger = sprintf("%.1f (%.1f - %.1f)", 
                                mean, 
                                min, 
                                max)) |>
  dplyr::select(Treatment, larger) |>
  tidyr::pivot_wider(
    names_from = Treatment,
    values_from = larger
    )
# A tibble: 1 × 5
  Full_Planted      Nuclei_25         Nuclei_50        Strip_25         Strip_50
  <chr>             <chr>             <chr>            <chr>            <chr>   
1 13.9 (0.0 - 35.0) 13.9 (0.0 - 39.0) 6.9 (0.0 - 43.0) 18.7 (0.0 - 104… 5.2 (0.…

Analysis for their abundances

show code
ruderal_df <- species_table |> 
  dplyr::filter(scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes"))|>
  # dplyr::select(Treatment, Block, Plot, Planted)|>unique()|>   with(table(Block, Plot, by=Treatment)) # Full_Planted P2B5 & P3B1 (Inside)
  dplyr::add_row(
    scientificName = c(NA, NA),
    `Seed dispersal` = c(NA, NA),
    `Planting group` = c(NA, NA),
    Treatment = c("Full_Planted", "Full_Planted"),
    Block = c("B5", "B1"),
    Plot = c("P2", "P3"),
    Planted = c("Inside","Inside"),
    n_sapling= c(0,0),
    n_tree = c(0,0),
    n_total= c(0,0)
    )|>
  dplyr::group_by(Treatment, Block, Plot, Planted)|>
  dplyr::summarise(n = sum(n_total), .groups="drop")

Plotting

show code
ruderal_df |>
  ggplot(aes(x=Planted, y=n, color=Block))+
  geom_point()+facet_wrap(~Treatment)+theme_bw()

Models

show code
modelo_rud1 <- glmmTMB::glmmTMB(n ~ Treatment*Planted + (1 | Block),
                   ziformula = ~1,
                    # family = poisson(),
                     family = glmmTMB::nbinom2(),
                    # family = glmmTMB::truncated_nbinom2(),
                   data = ruderal_df)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
show code
modelo_rud2 <-  glmmTMB::glmmTMB(n ~ Treatment + (1 | Block),
                   ziformula = ~1,
                   #  family = poisson(),
                    family = glmmTMB::nbinom2(),
                    # family = glmmTMB::truncated_nbinom2(),
                   data = ruderal_df)
modelo_rud3 <-  glmmTMB::glmmTMB(n ~ Planted + (1 | Block),
                   ziformula = ~1,
                   # family = poisson(),
                    family = glmmTMB::nbinom2(),
                  #   family = glmmTMB::truncated_nbinom2(),
                   data = ruderal_df)
modelo_rud_null <-  glmmTMB::glmmTMB(n ~ 0 + (1 | Block),
                    ziformula = ~1,
                    # family = poisson(), # ZIP
                     family = glmmTMB::nbinom2(), # ZINB
                    # family = glmmTMB::truncated_nbinom2(), # hurdle
                   data = ruderal_df)
# Model (null) | (treat*planted) | (treat) | (planted)
# ZIP

AIC(modelo_rud_null, modelo_rud1, modelo_rud2, modelo_rud3)
                df      AIC
modelo_rud_null  3 1725.166
modelo_rud1     12 1658.610
modelo_rud2      8 1707.148
modelo_rud3      5 1660.457
show code
performance::check_overdispersion(modelo_rud1)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
# Overdispersion test

 dispersion ratio = 0.965
          p-value = 0.872
No overdispersion detected.
show code
performance::check_model(modelo_rud1)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
`check_outliers()` does not yet support models of class `glmmTMB`.

Exploring residuals

show code
par(mfrow=c(2,2))
DHARMa::plotResiduals(modelo_rud1)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
show code
DHARMa::plotQQunif(modelo_rud1)
dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
show code
DHARMa::plotResiduals(modelo_rud3)
DHARMa::plotQQunif(modelo_rud3)

show code
car::Anova(modelo_rud1)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: n
                    Chisq Df Pr(>Chisq)    
Treatment         16.1644  4   0.002806 ** 
Planted           65.9720  1  4.574e-16 ***
Treatment:Planted  1.5434  3   0.672287    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Is there any difference between treatments considering inside/outside?

show code
emmeans::emmeans(modelo_rud1, ~ Treatment|Planted) |> 
  multcomp::cld(adjust = "bonferroni")
Planted = Inside:
 Treatment    emmean    SE  df asymp.LCL asymp.UCL .group
 Full_Planted   3.66 0.206 Inf      3.13      4.19  1    
 Strip_25       4.09 0.201 Inf      3.57      4.60  12   
 Nuclei_25      4.13 0.201 Inf      3.61      4.65  12   
 Nuclei_50      4.23 0.200 Inf      3.72      4.75   2   
 Strip_50       4.29 0.200 Inf      3.78      4.81   2   

Planted = Outside:
 Treatment    emmean    SE  df asymp.LCL asymp.UCL .group
 Strip_25       3.42 0.203 Inf      2.91      3.93  1    
 Strip_50       3.48 0.202 Inf      2.97      3.99  1    
 Nuclei_50      3.51 0.202 Inf      3.00      4.01  1    
 Nuclei_25      3.60 0.204 Inf      3.09      4.11  1    
 Full_Planted nonEst    NA  NA        NA        NA       

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for varying numbers of estimates 
P value adjustment: bonferroni method for varying numbers of tests 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Answer: Yes (planted) and no (only unplanted areas)

Is there any difference between inside/outside considering treatments?

show code
emmeans::emmeans(modelo_rud1, ~ Planted|Treatment) |> 
  multcomp::cld(adjust = "bonferroni")
Treatment = Full_Planted:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Inside    3.66 0.206 Inf      3.26      4.06  1    
 Outside nonEst    NA  NA        NA        NA       

Treatment = Nuclei_25:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Outside   3.60 0.204 Inf      3.15      4.06  1    
 Inside    4.13 0.201 Inf      3.68      4.58   2   

Treatment = Nuclei_50:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Outside   3.51 0.202 Inf      3.05      3.96  1    
 Inside    4.23 0.200 Inf      3.79      4.68   2   

Treatment = Strip_25:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Outside   3.42 0.203 Inf      2.97      3.88  1    
 Inside    4.09 0.201 Inf      3.64      4.54   2   

Treatment = Strip_50:
 Planted emmean    SE  df asymp.LCL asymp.UCL .group
 Outside   3.48 0.202 Inf      3.03      3.93  1    
 Inside    4.29 0.200 Inf      3.85      4.74   2   

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for varying numbers of estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Answer: No.

Table S3

show code
sjPlot::tab_model(modelo_rud1)
Model matrix is rank deficient. Parameters
  `TreatmentStrip_50:PlantedOutside` were not estimable.
  n
Predictors Incidence Rate Ratios CI p
Count Model
(Intercept) 38.84 25.95 – 58.15 <0.001
Treatment [Nuclei_25] 1.60 1.14 – 2.24 0.006
Treatment [Nuclei_50] 1.78 1.27 – 2.49 0.001
Treatment [Strip_25] 1.53 1.09 – 2.15 0.013
Treatment [Strip_50] 1.89 1.35 – 2.64 <0.001
Planted [Outside] 0.44 0.32 – 0.61 <0.001
Treatment [Nuclei_25] ×
Planted [Outside]
1.33 0.84 – 2.12 0.228
Treatment [Nuclei_50] ×
Planted [Outside]
1.09 0.69 – 1.73 0.711
Treatment [Strip_25] ×
Planted [Outside]
1.16 0.73 – 1.85 0.529
(Intercept) 51.38 23.07 – 140.38
Zero-Inflated Model
(Intercept) 0.01 0.00 – 0.05 <0.001
Random Effects
σ2 0.02
τ00 Block 0.13
ICC 0.87
N Block 5
Observations 179
Marginal R2 / Conditional R2 0.424 / 0.924
show code
broom.mixed::tidy(modelo_rud1) |>
  dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3))) |>
  dplyr::select(-c(component, group))|>
  knitr::kable()
effect term estimate std.error statistic p.value
fixed (Intercept) 3.659 0.206 17.777 0.000
fixed TreatmentNuclei_25 0.471 0.171 2.748 0.006
fixed TreatmentNuclei_50 0.575 0.171 3.356 0.001
fixed TreatmentStrip_25 0.427 0.172 2.486 0.013
fixed TreatmentStrip_50 0.635 0.171 3.713 0.000
fixed PlantedOutside -0.814 0.167 -4.873 0.000
fixed TreatmentNuclei_25:PlantedOutside 0.286 0.237 1.207 0.228
fixed TreatmentNuclei_50:PlantedOutside 0.087 0.236 0.371 0.711
fixed TreatmentStrip_25:PlantedOutside 0.149 0.237 0.630 0.529
fixed TreatmentStrip_50:PlantedOutside NA NA NA NA
fixed (Intercept) -4.497 0.721 -6.238 0.000
ran_pars sd__(Intercept) 0.365 NA NA NA

(iv) Figure Recruitment

Small letters: Comparison between treatments, only inside planted areas.

Capital letters: Comparison between treatments, only outside planted areas.

Smaller recruits

show code
ordem <- c("Nuclei_25", "Strip_25", "Nuclei_50", "Strip_50", "Full_Planted")

recruits_ha <- species_table |> 
  dplyr::mutate(
    n_small = dplyr::if_else(is.na(n_sapling), 0, n_sapling),
    n_large = dplyr::if_else(is.na(n_tree), 0, n_tree)) |>
  dplyr::group_by(Treatment, Block, Planted) |>
  dplyr::summarise(
    n_small_m2 = (sum(n_small, na.rm=T))/(144*4), 
    n_large_m2 = (sum(n_large, na.rm=T))/(144*4), 
    .groups = "drop") |>
  dplyr::mutate(
    n_small_ha = n_small_m2 * 10000,
    n_large_ha = n_large_m2 * 10000
  ) |>
  tidyr::pivot_wider(
    names_from = Planted,
    values_from = c(n_small_ha, n_large_ha),
    names_sep = "_"
  ) |>
  dplyr::select(Treatment, Block, n_small_ha_Inside:n_large_ha_Outside)


exc_recruits_ha <- species_table |> 
  dplyr::filter(!scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes")) |>
  dplyr::mutate(
    n_small = dplyr::if_else(is.na(n_sapling), 0, n_sapling),
    n_large = dplyr::if_else(is.na(n_tree), 0, n_tree)) |>
  dplyr::group_by(Treatment, Block, Planted) |>
  dplyr::summarise(
    n_small_m2 = (sum(n_small, na.rm=T))/(144*4), 
    n_large_m2 = (sum(n_large, na.rm=T))/(144*4), 
    .groups = "drop") |>
  dplyr::mutate(
    n_small_ha = n_small_m2 * 10000,
    n_large_ha = n_large_m2 * 10000
  ) |>
  tidyr::pivot_wider(
    names_from = Planted,
    values_from = c(n_small_ha, n_large_ha),
    names_sep = "_"
  ) |>
  dplyr::select(Treatment, Block, n_small_ha_Inside:n_large_ha_Outside)

ruderal_ha <- species_table |> 
  dplyr::filter(scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes")) |>
  dplyr::mutate(
    n_small = dplyr::if_else(is.na(n_sapling), 0, n_sapling),
    n_large = dplyr::if_else(is.na(n_tree), 0, n_tree)) |>
  dplyr::group_by(Treatment, Block, Planted) |>
  dplyr::summarise(
    n_small_m2 = (sum(n_small, na.rm=T))/(144*4), 
    n_large_m2 = (sum(n_large, na.rm=T))/(144*4), 
    .groups = "drop") |>
  dplyr::mutate(
    n_small_ha = n_small_m2 * 10000,
    n_large_ha = n_large_m2 * 10000
  ) |>
  tidyr::pivot_wider(
    names_from = Planted,
    values_from = c(n_small_ha, n_large_ha),
    names_sep = "_"
  ) |>
  dplyr::select(Treatment, Block, n_small_ha_Inside:n_large_ha_Outside)


prop_table <- tidyr::tibble(
  Treatment = unique(sort(species_table$Treatment)),
  prop_P = c(1, 0.25, 0.5, 0.25, 0.5),
  prop_UP = 1 - prop_P
)

data_fig_recruits <- recruits_ha |>
  dplyr::left_join(prop_table, by = "Treatment") |>
  dplyr::mutate(
    n_small_total_ha = dplyr::coalesce(n_small_ha_Inside * prop_P, 0)  + dplyr::coalesce(n_small_ha_Outside * prop_UP, 0),
    n_large_total_ha = dplyr::coalesce(n_large_ha_Inside * prop_P, 0) + dplyr::coalesce(n_large_ha_Outside * prop_UP, 0)
  ) |>
  dplyr::select(Treatment,
         starts_with("n_small_ha"),
         starts_with("n_large_ha"),
         n_small_total_ha,
         n_large_total_ha
  )
data_fig_ruderal <- ruderal_ha |>
  dplyr::left_join(prop_table, by = "Treatment") |>
  dplyr::mutate(
    n_small_total_ha = dplyr::coalesce(n_small_ha_Inside * prop_P, 0)  + dplyr::coalesce(n_small_ha_Outside * prop_UP, 0),
    n_large_total_ha = dplyr::coalesce(n_large_ha_Inside * prop_P, 0) + dplyr::coalesce(n_large_ha_Outside * prop_UP, 0)
  ) |>
  dplyr::select(Treatment,
         starts_with("n_small_ha"),
         starts_with("n_large_ha"),
         n_small_total_ha,
         n_large_total_ha
  )

data_fig_exc_ruderal <- exc_recruits_ha |>
  dplyr::left_join(prop_table, by = "Treatment") |>
  dplyr::mutate(
    n_small_total_ha = dplyr::coalesce(n_small_ha_Inside * prop_P, 0)  + dplyr::coalesce(n_small_ha_Outside * prop_UP, 0),
    n_large_total_ha = dplyr::coalesce(n_large_ha_Inside * prop_P, 0) + dplyr::coalesce(n_large_ha_Outside * prop_UP, 0)
  ) |>
  dplyr::select(Treatment,
         starts_with("n_small_ha"),
         starts_with("n_large_ha"),
         n_small_total_ha,
         n_large_total_ha
  )

Small regenerants

show code
data_fig_small_1 <- data_fig_recruits |>
  dplyr::group_by(Treatment)|>
  dplyr::summarise(
    mean_inside = round(mean(n_small_ha_Inside, na.rm=TRUE), 1),
    se_inside = round((sd(n_small_ha_Inside, na.rm=TRUE))/sqrt(dplyr::n()), 1),
    mean_outside = round(mean(n_small_ha_Outside, na.rm=TRUE), 1),
    se_outside = round((sd(n_small_ha_Outside, na.rm=TRUE))/sqrt(dplyr::n()), 1)
  ) |>
   tidyr::pivot_longer(cols = -Treatment, names_to = c(".value", "side"), names_sep = "_") |> 
  dplyr::filter(!is.na(se)) |>
  dplyr::mutate(class="all recruits")

data_fig_small_2 <- data_fig_ruderal |>
  dplyr::group_by(Treatment)|>
  dplyr::summarise(
    mean_inside = round(mean(n_small_ha_Inside, na.rm=TRUE), 1),
    se_inside = round((sd(n_small_ha_Inside, na.rm=TRUE))/sqrt(dplyr::n()), 1),
    mean_outside = round(mean(n_small_ha_Outside, na.rm=TRUE), 1),
    se_outside = round((sd(n_small_ha_Outside, na.rm=TRUE))/sqrt(dplyr::n()), 1)
  ) |>
   tidyr::pivot_longer(cols = -Treatment, names_to = c(".value", "side"), names_sep = "_") |> 
  dplyr::filter(!is.na(se)) |>
  dplyr::mutate(class="ruderal")

data_fig_small_3 <- data_fig_exc_ruderal |>
  dplyr::group_by(Treatment)|>
  dplyr::summarise(
    mean_inside = round(mean(n_small_ha_Inside, na.rm=TRUE), 1),
    se_inside = round((sd(n_small_ha_Inside, na.rm=TRUE))/sqrt(dplyr::n()), 1),
    mean_outside = round(mean(n_small_ha_Outside, na.rm=TRUE), 1),
    se_outside = round((sd(n_small_ha_Outside, na.rm=TRUE))/sqrt(dplyr::n()), 1)
  ) |>
   tidyr::pivot_longer(cols = -Treatment, names_to = c(".value", "side"), names_sep = "_") |> 
  dplyr::filter(!is.na(se)) |>
  dplyr::mutate(class = "except ruderal")

data_fig_small_alterado <- data_fig_small_1

data_fig_small_alterado$mean <- data_fig_small_3$mean

data_fig_small <- dplyr::bind_rows(data_fig_small_alterado, data_fig_small_2)

New Figure

show code
ordem_x <- c(
  "Nuclei_25_inside", "Nuclei_25_outside",
  "Strip_25_inside",  "Strip_25_outside",
  "Nuclei_50_inside", "Nuclei_50_outside",
  "Strip_50_inside",  "Strip_50_outside",
  "Full_Planted_inside"
)

data_fig_small_top2 <- data_fig_small |>
  group_by(Treatment, side) |>
  arrange(desc(class)) |>
  mutate(cum_mean = cumsum(mean)) |>
  ungroup() |>
  mutate(Treat_side = interaction(Treatment, side, sep = "_"),
         Treat_side = factor(Treat_side, levels = ordem_x))|>
  dplyr::arrange(desc(class),Treat_side)

letras <- data.frame(
  Treat_side = levels(data_fig_small_top2$Treat_side),
  label = c("AB","a",
            "AB","a",
            "B","a",
            "B","a",
            "A"),  
  y = c(5500, 5000, 5200, 4500, 
        7000, 4500, 7000, 4800, 3800)           
)


fig_small<-ggplot(data_fig_small_top2, 
       aes(x = Treat_side, y = mean, fill = side, pattern = class)) +
  geom_col_pattern(
    position = "stack",
    color = "black",
    pattern_fill = "#92D050", pattern_alpha = 0.5,
    pattern_angle = 45, pattern_density = 0.1,
    pattern_spacing = 0.1
  ) +
  scale_pattern_manual(
    values = c("regenerant" = "none", "ruderal" = "stripe"),
    guide = "none") +
  scale_fill_manual(labels = c("Planted", "Unplanted"), 
                    values = c("#92D050", "white")) +
  geom_errorbar(
    data = data_fig_small_top2 |> group_by(Treatment, side) |> slice_tail(n = 1),
    aes(ymin = cum_mean - se, ymax = cum_mean + se, x = Treat_side),
    width = .2,
    inherit.aes = FALSE
  ) +
  labs(fill = "", 
       y =  expression("Recruits (n.ha"^{-1}~")"),
       x = "") +
  geom_text(data = letras, 
          aes(x = Treat_side, y = y, label = label), 
          inherit.aes = FALSE,  
          size = 3.5, vjust = -0.5) +
  scale_x_discrete(labels = function(x) gsub("_", "\n", x)) +
  scale_y_continuous(expand =c(0,0), limits = c(0,8000), breaks = seq(0,  8000, by=1500))+
  theme_classic() +
  theme(
    axis.text.x = element_blank(), #element_text(color="black", hjust = 0.5, size=9),
    axis.text.y = element_text(color="black", size=10),
    axis.title.y = element_text(color="black", size=11),
    legend.position = c(0.18, 0.95),
    legend.text = element_text(size = 10),
    legend.key.size = unit(3, "mm"),
    legend.spacing = unit(0, "mm"),
    legend.background = element_blank(),
    legend.title = element_text(size = 8))
fig_small

show code
# ggsave("fig_small_recruits.png", plot=fig_small, width=90, height = 90 ,unit="mm")

Larger recruits

show code
data_fig_large <- result_regenerant_all |>
  dplyr::group_by(Treatment)|>
  dplyr::summarise(
    mean_inside = round(mean(n_large_ha_Inside, na.rm=TRUE), 1),
    se_inside = round((sd(n_large_ha_Inside, na.rm=TRUE))/sqrt(dplyr::n()), 1),
    mean_outside = round(mean(n_large_ha_Outside, na.rm=TRUE), 1),
    se_outside = round((sd(n_large_ha_Outside, na.rm=TRUE))/sqrt(dplyr::n()), 1)
  ) |>
   tidyr::pivot_longer(cols = -Treatment, names_to = c(".value", "side"), names_sep = "_") |> 
  dplyr::mutate(
    Treatment = ordered(Treatment,levels = ordem),
    tratamento = forcats::fct_recode(Treatment, 
                                     "Nuclei 25%" = "Nuclei_25", "Strips 25%" = "Strip_25", 
                                     "Nuclei 50%" = "Nuclei_50", "Strips 50%" = "Strip_50","Planting 100%" = "Full_Planted"),
    Treat_side = interaction(Treatment, side, sep = "_"),
    Treat_side = factor(Treat_side, ordered = TRUE, levels = ordem_x,
                         labels = c("Nuclei\n  25%", " " ,
                                    "Strips\n  25%", "  ", 
                                    "Nuclei\n  50%", "   ",
                                    "Strips\n  50%", "    ",
                                    "Planting\n100%"))
  )|>
  dplyr::filter(!is.na(se))

Another way

show code
letras <- data.frame(
  Treat_side = levels(data_fig_large$Treat_side),
  label = c("A","a",
            "A","a",
            "A","a",
            "A","a",
            "A"),  
  y = c(310, 470, 
        200, 370, 
        200, 200, 
        300, 470, 
        320)            
)


data_fig_large <- data_fig_large |>
  dplyr::mutate(hjust_value = ifelse(Treat_side == last(Treat_side), 0.5, 0))

fig_large<-ggplot(data_fig_large, 
       aes(x = Treat_side, y = mean, fill = side)) +
  geom_bar(
    stat = "identity",
    position = "stack",
    color = "black") +
  scale_fill_manual(labels = c("Planted", "Unplanted"), 
                    values = c("#92D050", "white")) +
  geom_errorbar(
    aes(ymin = mean - se, ymax = mean + se, x = Treat_side),
    width = .2
  ) +
  labs(fill = "", 
       y =  expression("Recruits (n.ha"^{-1}*")"),
       x = "") +
  geom_text(data = letras, aes(x = Treat_side, y = y, label = label),  inherit.aes = FALSE,  size = 3.5, vjust = -0.5) +
  scale_y_continuous(expand =c(0,0), limits = c(0,550), breaks = seq(0, 500, by=100))+
  theme_classic() +
  theme(
    axis.text.x = element_text(color="black", hjust = data_fig_large$hjust_value, size=9),
    axis.text.y = element_text(color="black", size=10),
    axis.title.y = element_text(color="black", size=11),
    legend.position = "none",# c(0.18, 0.95),
    legend.text = element_text(size = 10),
    legend.key.size = unit(3, "mm"),
    legend.spacing = unit(0, "mm"),
    legend.background = element_blank(),
    legend.title = element_text(size = 8))
Warning: Vectorized input to `element_text()` is not officially supported.
ℹ Results may be unexpected or may change in future versions of ggplot2.
show code
fig_large

Altogether:

show code
fig_recruits <- cowplot::plot_grid(fig_small + theme(axis.text.x=element_blank()), 
                                   fig_large, 
                                   ncol = 1,
                                   labels = c("(a)", "(b)"), label_size = 10, label_x = 0.01, label_y = c(1,1.1)
                                   )
fig_recruits

show code
ggsave("fig_recruits.png", plot=fig_recruits, width=90, height = 100 ,unit="mm")

(vi) Species list table

Only recruits (not including planted individuals)

show code
table_recruits <- readxl::read_xlsx("Itatinga_Nuc_recruitment_SpeciesList.xlsx") 
table_recruits |> knitr::kable(format = "html")
family scientificName Nuclei 25% (small) Nuclei 25% (large) Nuclei 25% (total) Strips 25% (small) Strips 25% (large) Strips 25% (total) Nuclei 50% (small) Nuclei 50% (large) Nuclei 50% (total) Strips 50% (small) Strips 50% (large) Strips 50% (total) Planting 100% (small) Planting 100% (large) Planting 100% (total) Total small Total large Total
Anacardiaceae Astronium graveolens 0 0 0 0 0 0 1 0 1 0 5 5 0 0 0 1 5 6
Anacardiaceae Myracrodruon urundeuva 0 12 12 0 2 2 0 2 2 0 0 0 0 0 0 0 16 16
Anacardiaceae Schinus terebinthifolia 0 0 0 1 0 1 1 0 1 0 0 0 0 0 0 2 0 2
Anacardiaceae Tapirira guianensis 0 3 3 0 0 0 0 0 0 0 0 0 0 0 0 0 3 3
Annonaceae Annona cacans 0 0 0 11 0 11 0 0 0 0 2 2 1 0 1 12 2 14
Annonaceae Annona coriacea 0 0 0 2 0 2 0 2 2 0 0 0 0 0 0 2 2 4
Annonaceae Duguetia furfuracea 0 0 0 1 0 1 0 0 0 0 1 1 1 0 1 2 1 3
Apocynaceae Aspidosperma tomentosum 2 1 3 0 0 0 0 0 0 0 0 0 0 4 4 2 5 7
Araliaceae Didymopanax vinosus 0 0 0 4 3 7 2 0 2 3 1 4 0 0 0 9 4 13
Arecaceae Syagrus romanzoffiana 1 1 2 1 0 1 3 0 3 0 0 0 0 0 0 5 1 6
Asteraceae Baccharis dracunculifolia 1432 14 1446 916 15 931 931 8 939 1784 6 1790 511 4 515 5574 47 5621
Asteraceae Moquiniastrum barrosoae 0 0 0 3 1 4 0 0 0 0 0 0 0 0 0 3 1 4
Asteraceae Moquiniastrum polymorphum 0 17 17 2 3 5 13 2 15 2 0 2 0 0 0 17 22 39
Asteraceae Piptocarpha rotundifolia 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1 0 1
Asteraceae Vernonanthura polyanthes 634 0 634 910 0 910 1311 0 1311 491 0 491 0 2 2 3346 2 3348
Bignoniaceae Cybistax antisyphilitica 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Bignoniaceae Handroanthus albus 0 2 2 0 0 0 0 0 0 0 0 0 9 0 9 9 2 11
Bignoniaceae Handroanthus ochraceus 1 0 1 0 0 0 0 2 2 0 0 0 0 0 0 1 2 3
Boraginaceae Cordia sellowiana 0 0 0 0 6 6 0 3 3 0 3 3 0 0 0 0 12 12
Burseraceae Protium spruceanum 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 1 1
Calophyllaceae Myrsine gardneriana 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1
Cannabaceae Trema micrantha 0 0 0 0 1 1 0 0 0 0 10 10 0 0 0 0 11 11
Chrysobalanaceae Couepia grandiflora 1 0 1 0 2 2 0 0 0 0 0 0 0 0 0 1 2 3
Combretaceae Terminalia argentea 0 0 0 0 2 2 0 2 2 0 0 0 0 0 0 0 4 4
Euphorbiaceae Alchornea triplinervia 0 0 0 0 0 0 0 0 0 1 0 1 44 0 44 45 0 45
Euphorbiaceae Croton floribundus 73 0 73 55 1 56 119 0 119 147 0 147 0 0 0 394 1 395
Euphorbiaceae Croton urucurana 1 0 1 2 0 2 0 0 0 7 1 8 0 0 0 10 1 11
Euphorbiaceae Mabea fistulifera 1 0 1 12 0 12 7 0 7 31 2 33 0 3 3 51 5 56
Euphorbiaceae Sapium glandulosum 1 0 1 0 0 0 0 0 0 0 2 2 0 0 0 1 2 3
Fabaceae Anadenanthera colubrina var. cebil 2 5 7 5 0 5 0 1 1 8 13 21 0 0 0 15 19 34
Fabaceae Andira fraxinifolia 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1 0 1
Fabaceae Bauhinia rufa 0 0 0 0 0 0 0 0 0 2 0 2 0 0 0 2 0 2
Fabaceae Chamaecrista comosa 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Fabaceae Copaifera langsdorffii 2 2 4 2 4 6 2 2 4 5 0 5 0 0 0 11 8 19
Fabaceae Dalbergia frutescens 0 2 2 7 0 7 0 0 0 0 0 0 0 0 0 7 2 9
Fabaceae Dalbergia miscolobium 0 0 0 1 2 3 1 1 2 2 1 3 0 0 0 4 4 8
Fabaceae Enterolobium contortisiliquum 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 1
Fabaceae Erythrina mulungu 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 1 0 1
Fabaceae Hymenaea stigonocarpa 0 2 2 0 1 1 0 0 0 0 0 0 0 0 0 0 3 3
Fabaceae Inga vera 0 0 0 4 0 4 4 0 4 3 0 3 1 0 1 12 0 12
Fabaceae Leptolobium dasycarpum 0 0 0 0 0 0 0 3 3 0 0 0 0 0 0 0 3 3
Fabaceae Leptolobium elegans 2 32 34 1 17 18 3 32 35 3 0 3 1 0 1 10 81 91
Fabaceae Leucochloron incuriale 0 0 0 0 5 5 0 0 0 0 2 2 0 0 0 0 7 7
Fabaceae Machaerium acutifolium 1 0 1 1 3 4 0 0 0 0 0 0 0 0 0 2 3 5
Fabaceae Machaerium brasiliense 1 0 1 1 5 6 0 0 0 1 0 1 0 0 0 3 5 8
Fabaceae Machaerium villosum 0 0 0 0 3 3 0 0 0 0 0 0 0 0 0 0 3 3
Fabaceae Platypodium elegans 0 0 0 0 0 0 1 0 1 0 0 0 1 0 1 2 0 2
Fabaceae Pterogyne nitens 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1
Fabaceae Senegalia polyphylla 48 0 48 94 1 95 87 2 89 101 0 101 0 0 0 330 3 333
Fabaceae Senna pendula 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 1 1
Fabaceae Senna rugosa 6 0 6 7 0 7 9 0 9 7 35 42 0 28 28 29 63 92
Lamiaceae Aegiphila verticillata 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1 0 1
Lamiaceae Vitex polygama 0 0 0 0 0 0 0 0 0 0 2 2 14 0 14 14 2 16
Lauraceae Aiouea sellowiana 0 2 2 0 0 0 0 2 2 0 0 0 0 0 0 0 4 4
Lauraceae Endlicheria paniculata 0 0 0 0 0 0 0 0 0 2 2 4 0 0 0 2 2 4
Lauraceae Ocotea corymbosa 1 0 1 1 0 1 0 1 1 0 0 0 0 0 0 2 1 3
Lauraceae Ocotea puberula 0 0 0 0 0 0 2 0 2 0 0 0 0 0 0 2 0 2
Lauraceae Ocotea pulchella 9 5 14 2 0 2 2 1 3 0 6 6 0 0 0 13 12 25
Lauraceae Ocotea velutina 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 1 1
Lythraceae Lafoensia pacari 0 0 0 1 1 2 0 0 0 0 0 0 0 0 0 1 1 2
Malpighiaceae Byrsonima crassifolia 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1
Malpighiaceae Byrsonima intermedia 0 0 0 6 1 7 3 0 3 0 6 6 0 0 0 9 7 16
Malvaceae Eriotheca gracilipes 1 1 2 0 0 0 0 0 0 0 0 0 0 0 0 1 1 2
Melastomataceae Miconia albicans 13 12 25 2 6 8 0 3 3 9 4 13 7 0 7 31 25 56
Melastomataceae Miconia auricoma 1 0 1 0 0 0 0 0 0 0 11 11 1 3 4 2 14 16
Melastomataceae Miconia langsdorffii 0 5 5 0 0 0 0 0 0 0 0 0 2 0 2 2 5 7
Melastomataceae Miconia ligustroides 5 0 5 5 0 5 0 0 0 6 7 13 0 2 2 16 9 25
Melastomataceae Miconia minutiflora 0 0 0 0 0 0 0 0 0 0 3 3 2 0 2 2 3 5
Myrtaceae Campomanesia pubescens 0 0 0 0 0 0 0 0 0 1 1 2 0 0 0 1 1 2
Myrtaceae Eugenia aurata 0 0 0 1 0 1 0 0 0 1 6 7 0 0 0 2 6 8
Myrtaceae Eugenia bimarginata 1 0 1 1 2 3 0 0 0 1 0 1 0 0 0 3 2 5
Myrtaceae Eugenia hyemalis 1 0 1 9 0 9 0 0 0 0 0 0 0 0 0 10 0 10
Myrtaceae Eugenia myrcianthes 0 0 0 2 1 3 0 0 0 0 0 0 0 0 0 2 1 3
Myrtaceae Eugenia punicifolia 5 0 5 6 0 6 0 0 0 1 0 1 0 0 0 12 0 12
Myrtaceae Eugenia pyriformis 1 0 1 6 0 6 0 0 0 6 0 6 0 0 0 13 0 13
Myrtaceae Eugenia uniflora 0 0 0 5 0 5 0 0 0 0 0 0 0 0 0 5 0 5
Myrtaceae Myrcia 0 0 0 0 0 0 2 0 2 0 0 0 0 0 0 2 0 2
Myrtaceae Myrcia bella 2 14 16 3 11 14 5 8 13 2 3 5 0 0 0 12 36 48
Myrtaceae Myrcia guianensis 5 5 10 1 0 1 4 2 6 0 0 0 0 0 0 10 7 17
Myrtaceae Myrcia rufipes 1 4 5 2 0 2 1 0 1 3 0 3 0 0 0 7 4 11
Myrtaceae Myrcia splendens 4 0 4 4 0 4 0 1 1 4 3 7 0 0 0 12 4 16
Myrtaceae Myrcia tomentosa 0 0 0 0 2 2 0 0 0 0 1 1 0 0 0 0 3 3
Myrtaceae Myrciaria floribunda 1 0 1 0 0 0 0 0 0 0 6 6 0 0 0 1 6 7
Myrtaceae Myrciaria glazioviana 11 0 11 4 1 5 5 0 5 7 0 7 0 0 0 27 1 28
Myrtaceae Psidium grandifolium 0 0 0 0 0 0 0 0 0 1 0 1 1 0 1 2 0 2
Myrtaceae Psidium guineense 1 0 1 1 0 1 0 0 0 0 0 0 0 0 0 2 0 2
Myrtaceae Psidium oligospermum 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1
Nyctaginaceae Guapira noxia 0 4 4 0 1 1 0 0 0 3 0 3 0 0 0 3 5 8
Ochnaceae Ouratea spectabilis 9 12 21 4 9 13 1 0 1 6 0 6 0 0 0 20 21 41
Peraceae Pera glabrata 0 0 0 1 0 1 4 8 12 1 1 2 0 2 2 6 11 17
Piperaceae Piper mollicomum 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1 0 1
Polygalaceae Bredemeyera floribunda 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1 0 1
Polygonaceae Coccoloba cordata 0 0 0 3 5 8 1 0 1 1 0 1 0 0 0 5 5 10
Polygonaceae Triplaris americana 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1
Primulaceae Myrsine coriacea 0 0 0 0 0 0 0 0 0 2 6 8 0 0 0 2 6 8
Primulaceae Myrsine gardneriana 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1
Primulaceae Myrsine guianensis 0 0 0 4 0 4 0 3 3 0 0 0 0 0 0 4 3 7
Primulaceae Myrsine umbellata 2 0 2 0 0 0 0 0 0 3 0 3 0 0 0 5 0 5
Proteaceae Roupala montana 1 2 3 0 3 3 0 0 0 3 0 3 0 0 0 4 5 9
Rubiaceae Amaioua intermedia 0 0 0 1 0 1 0 0 0 4 0 4 18 0 18 23 0 23
Rubiaceae Palicourea sessilis 1 0 1 0 0 0 0 0 0 4 1 5 0 0 0 5 1 6
Salicaceae Casearia sylvestris 36 0 36 76 0 76 60 0 60 15 0 15 2 0 2 189 0 189
Sapindaceae Matayba elaeagnoides 0 0 0 3 0 3 2 0 2 0 0 0 0 0 0 5 0 5
Sapindaceae Sapindus saponaria 0 0 0 1 0 1 0 0 0 0 3 3 11 0 11 12 3 15
Sapotaceae Pouteria ramiflora 0 0 0 0 0 0 0 1 1 0 1 1 26 0 26 26 2 28
Solanaceae Cestrum corymbosum 0 0 0 0 0 0 0 0 0 0 3 3 3 0 3 3 3 6
Solanaceae Solanum didymum 0 0 0 2 0 2 0 0 0 0 0 0 0 0 0 2 0 2
Solanaceae Solanum granuloso-leprosum 44 3 47 79 6 85 77 1 78 55 3 58 2 7 9 257 20 277
Solanaceae Solanum sisymbriifolium 166 0 166 72 0 72 80 0 80 98 0 98 0 0 0 416 0 416
Urticaceae Cecropia pachystachya 1 0 1 2 0 2 1 0 1 0 0 0 0 9 9 4 9 13
Verbenaceae Lantana camara 0 0 0 2 0 2 0 0 0 0 12 12 1 0 1 3 12 15
Vochysiaceae Qualea grandiflora 0 0 0 0 0 0 2 0 2 0 0 0 222 0 222 224 0 224
NA Unclassified 9 8 17 5 4 9 8 1 9 12 2 14 0 0 0 34 15 49

! Altogether

show code
sampled_species <-readxl::read_xlsx("itatinga_sampled_species.xlsx")
sampled_species |> 
  dplyr::select(-`Seed dispersal`) |>
  knitr::kable(format = "html")
family scientificName scientificNameAuthorship Large_recruit Small_recruit Origin Planting group Successional group
Anacardiaceae Astronium graveolens Jacq. TRUE TRUE planted diversity non-pioneer
Anacardiaceae Schinus terebinthifolia Raddi FALSE TRUE recruit NA pioneer
Anacardiaceae Myracrodruon urundeuva Allem. TRUE FALSE recruit NA non-pioneer
Anacardiaceae Tapirira guianensis Aubl. TRUE FALSE recruit NA non-pioneer
Annonaceae Annona cacans Warm. TRUE TRUE planted diversity non-pioneer
Annonaceae Annona coriacea Mart. TRUE TRUE recruit NA non-pioneer
Annonaceae Duguetia furfuracea (A.St.-Hil.) Saff. TRUE TRUE recruit NA non-pioneer
Apocynaceae Aspidosperma tomentosum Mart. TRUE TRUE recruit NA non-pioneer
Araliaceae Didymopanax vinosus (Cham. & Schltdl.) Marchal TRUE TRUE recruit NA non-pioneer
Arecaceae Syagrus romanzoffiana (Cham.) Glassman TRUE TRUE recruit NA non-pioneer
Asteraceae Piptocarpha rotundifolia Baker FALSE TRUE recruit NA pioneer
Asteraceae Baccharis dracunculifolia DC. TRUE TRUE recruit NA pioneer
Asteraceae Moquiniastrum barrosoae (Cabrera) G.Sancho TRUE TRUE recruit NA pioneer
Asteraceae Moquiniastrum polymorphum (Less.) G.Sancho TRUE TRUE recruit NA pioneer
Asteraceae Vernonanthura polyanthes (Spreng.) A.J.Vega & Dematt. TRUE TRUE recruit NA pioneer
Bignoniaceae Handroanthus ochraceus (Cham.) Mattos TRUE TRUE planted diversity non-pioneer
Bignoniaceae Handroanthus impetiginosus (Mart. ex DC.) Mattos FALSE FALSE planted diversity non-pioneer
Bignoniaceae Jacaranda cuspidifolia Mart. FALSE FALSE planted diversity non-pioneer
Bignoniaceae Cybistax antisyphilitica (Mart.) Mart. TRUE FALSE recruit NA non-pioneer
Bignoniaceae Handroanthus albus (Cham.) Mattos TRUE TRUE recruit NA non-pioneer
Boraginaceae Cordia ecalyculata Vell. FALSE FALSE planted diversity pioneer
Boraginaceae Cordia sellowiana Cham. TRUE FALSE recruit NA pioneer
Burseraceae Protium spruceanum Engl. TRUE FALSE recruit NA non-pioneer
Cannabaceae Trema micrantha (L.) Blume TRUE FALSE recruit NA pioneer
Chrysobalanaceae Couepia grandiflora (Mart. & Zucc.) Benth. ex Hook.f. TRUE TRUE recruit NA non-pioneer
Combretaceae Terminalia argentea Mart. TRUE FALSE recruit NA non-pioneer
Euphorbiaceae Croton floribundus Spreng. TRUE TRUE planted filling pioneer
Euphorbiaceae Croton urucurana Baill. TRUE TRUE planted filling pioneer
Euphorbiaceae Mabea fistulifera Mart. TRUE TRUE planted filling pioneer
Euphorbiaceae Alchornea triplinervia (Spreng.) Müll.Arg. FALSE TRUE recruit NA non-pioneer
Euphorbiaceae Sapium glandulosum (L.) Morong TRUE TRUE recruit NA non-pioneer
Fabaceae Enterolobium contortisiliquum (Vell.) Morong FALSE TRUE planted diversity non-pioneer
Fabaceae Erythrina mulungu Mart. ex Benth. FALSE TRUE planted diversity non-pioneer
Fabaceae Platypodium elegans Vogel FALSE TRUE planted diversity non-pioneer
Fabaceae Copaifera langsdorffii Desf. TRUE TRUE planted diversity non-pioneer
Fabaceae Machaerium acutifolium Vogel TRUE TRUE planted diversity non-pioneer
Fabaceae Pterogyne nitens Tul. TRUE FALSE planted diversity non-pioneer
Fabaceae Anadenanthera colubrina var. cebil (Griseb.) Altschul FALSE FALSE planted diversity non-pioneer
Fabaceae Anadenanthera peregrina (L.) Speg. FALSE FALSE planted diversity non-pioneer
Fabaceae Apuleia leiocarpa (Vogel) J.F.Macbr. FALSE FALSE planted diversity non-pioneer
Fabaceae Cassia fistula L. FALSE FALSE planted diversity non-pioneer
Fabaceae Hymenaea courbaril L. FALSE FALSE planted diversity non-pioneer
Fabaceae Myroxylon peruiferum L.f. FALSE FALSE planted diversity non-pioneer
Fabaceae Parapiptadenia rigida (Benth.) Brenan FALSE FALSE planted diversity non-pioneer
Fabaceae Peltophorum dubium (Spreng.) Taub. FALSE FALSE planted diversity non-pioneer
Fabaceae Senna macranthera (DC. ex Collad.) H.S.Irwin & Barneby FALSE FALSE planted diversity non-pioneer
Fabaceae Inga vera Willd. FALSE TRUE planted filling non-pioneer
Fabaceae Senegalia polyphylla (DC.) Britton & Rose TRUE TRUE planted filling pioneer
Fabaceae Inga edulis Mart. FALSE FALSE planted filling non-pioneer
Fabaceae Andira fraxinifolia Benth. FALSE TRUE recruit NA non-pioneer
Fabaceae Bauhinia rufa (Bong.) Steud. FALSE TRUE recruit NA non-pioneer
Fabaceae Chamaecrista comosa E.Mey. TRUE FALSE recruit NA non-pioneer
Fabaceae Dalbergia frutescens (Vell.) Britton TRUE TRUE recruit NA non-pioneer
Fabaceae Dalbergia miscolobium Benth. TRUE TRUE recruit NA non-pioneer
Fabaceae Hymenaea stigonocarpa Mart. ex Hayne TRUE FALSE recruit NA non-pioneer
Fabaceae Leptolobium dasycarpum Vogel TRUE FALSE recruit NA non-pioneer
Fabaceae Leptolobium elegans Vogel TRUE TRUE recruit NA non-pioneer
Fabaceae Leucochloron incuriale (Vell.) Barneby & J.W.Grimes TRUE FALSE recruit NA non-pioneer
Fabaceae Machaerium brasiliense Vogel TRUE TRUE recruit NA non-pioneer
Fabaceae Machaerium villosum Vogel TRUE FALSE recruit NA non-pioneer
Fabaceae Senna pendula (Humb. & Bonpl. ex Willd.) H.S.Irwin & Barneby TRUE FALSE recruit NA pioneer
Fabaceae Senna rugosa (G.Don) H.S.Irwin & Barneby TRUE TRUE recruit NA pioneer
Lamiaceae Aegiphila verticillata Vell. FALSE TRUE recruit NA pioneer
Lamiaceae Vitex polygama Cham. TRUE TRUE recruit NA non-pioneer
Lauraceae Ocotea puberula (Rich.) Nees FALSE TRUE recruit NA non-pioneer
Lauraceae Aiouea sellowiana (Nees & Mart.) R.Rohde TRUE FALSE recruit NA non-pioneer
Lauraceae Endlicheria paniculata (Spreng.) J.F.Macbr. TRUE TRUE recruit NA non-pioneer
Lauraceae Ocotea corymbosa (Meisn.) Mez TRUE TRUE recruit NA non-pioneer
Lauraceae Ocotea pulchella (Nees & Mart.) Mez TRUE TRUE recruit NA non-pioneer
Lauraceae Ocotea velutina (Nees) Mart. ex B.D.Jacks. TRUE FALSE recruit NA non-pioneer
Lecythidaceae Cariniana estrellensis (Raddi) Kuntze FALSE FALSE planted diversity non-pioneer
Lythraceae Lafoensia pacari A.St.-Hil. TRUE TRUE planted diversity non-pioneer
Malpighiaceae Byrsonima crassifolia Kunth FALSE TRUE recruit NA non-pioneer
Malpighiaceae Byrsonima intermedia A.Juss. TRUE TRUE recruit NA non-pioneer
Malvaceae Ceiba speciosa (A.St.-Hil., A.Juss. & Cambess.) Ravenna FALSE FALSE planted diversity non-pioneer
Malvaceae Heliocarpus popayanensis Kunth FALSE FALSE planted diversity pioneer
Malvaceae Luehea divaricata Mart. FALSE FALSE planted diversity non-pioneer
Malvaceae Guazuma ulmifolia Lam. FALSE FALSE planted filling pioneer
Malvaceae Eriotheca gracilipes (K.Schum.) A.Robyns TRUE TRUE recruit NA non-pioneer
Melastomataceae Miconia albicans (Sw.) Steud. TRUE TRUE recruit NA pioneer
Melastomataceae Miconia auricoma (Spring ex Mart.) R.Goldenb. TRUE TRUE recruit NA pioneer
Melastomataceae Miconia langsdorffii Cogn. TRUE TRUE recruit NA pioneer
Melastomataceae Miconia ligustroides Naudin TRUE TRUE recruit NA pioneer
Melastomataceae Miconia minutiflora (Bonpl.) DC. TRUE TRUE recruit NA pioneer
Meliaceae Cedrela fissilis Vell. FALSE FALSE planted diversity non-pioneer
Moraceae Ficus guaranitica Chodat FALSE FALSE planted diversity non-pioneer
Myrtaceae Eugenia uniflora L. FALSE TRUE planted diversity non-pioneer
Myrtaceae Psidium myrtoides O.Berg FALSE FALSE planted diversity non-pioneer
Myrtaceae Eugenia hyemalis Cambess. FALSE TRUE recruit NA non-pioneer
Myrtaceae Eugenia punicifolia (Kunth) DC. FALSE TRUE recruit NA non-pioneer
Myrtaceae Eugenia pyriformis Cambess. FALSE TRUE recruit NA non-pioneer
Myrtaceae Myrcia DC. FALSE TRUE recruit NA non-pioneer
Myrtaceae Psidium grandifolium Mart. ex DC. FALSE TRUE recruit NA non-pioneer
Myrtaceae Psidium guineense Sw. FALSE TRUE recruit NA non-pioneer
Myrtaceae Campomanesia pubescens O.Berg TRUE TRUE recruit NA non-pioneer
Myrtaceae Eugenia aurata O.Berg TRUE TRUE recruit NA non-pioneer
Myrtaceae Eugenia bimarginata DC. TRUE TRUE recruit NA non-pioneer
Myrtaceae Eugenia myrcianthes Nied. TRUE TRUE recruit NA non-pioneer
Myrtaceae Myrcia bella Cambess. TRUE TRUE recruit NA non-pioneer
Myrtaceae Myrcia guianensis DC. TRUE TRUE recruit NA non-pioneer
Myrtaceae Myrcia rufipes DC. TRUE TRUE recruit NA non-pioneer
Myrtaceae Myrcia splendens DC. TRUE TRUE recruit NA non-pioneer
Myrtaceae Myrcia tomentosa DC. TRUE FALSE recruit NA non-pioneer
Myrtaceae Myrciaria floribunda O.Berg TRUE TRUE recruit NA non-pioneer
Myrtaceae Myrciaria glazioviana (Kiaersk.) G.M.Barroso ex Sobral TRUE TRUE recruit NA non-pioneer
Myrtaceae Psidium oligospermum Mart. ex DC. TRUE FALSE recruit NA non-pioneer
Nyctaginaceae Guapira noxia (Netto) Lundell TRUE TRUE recruit NA non-pioneer
Ochnaceae Ouratea spectabilis Engl. TRUE TRUE recruit NA non-pioneer
Peraceae Pera glabrata (Schott) Poepp. ex Baill. TRUE TRUE recruit NA non-pioneer
Petiveriaceae Gallesia integrifolia (Spreng.) Harms FALSE FALSE planted diversity non-pioneer
Piperaceae Piper mollicomum Kunth FALSE TRUE recruit NA non-pioneer
Polygalaceae Bredemeyera floribunda Willd. FALSE TRUE recruit NA non-pioneer
Polygonaceae Triplaris americana L. FALSE TRUE recruit NA non-pioneer
Polygonaceae Coccoloba cordata Cham. TRUE TRUE recruit NA non-pioneer
Primulaceae Myrsine gardneriana A.DC. FALSE TRUE recruit NA non-pioneer
Primulaceae Myrsine gardneriana A.DC. FALSE TRUE recruit NA non-pioneer
Primulaceae Myrsine umbellata Mart. FALSE TRUE recruit NA non-pioneer
Primulaceae Myrsine coriacea (Sw.) R.Br. ex Roem. & Schult. TRUE TRUE recruit NA non-pioneer
Primulaceae Myrsine guianensis (Aubl.) Kuntze TRUE TRUE recruit NA non-pioneer
Proteaceae Roupala montana Aubl. TRUE TRUE recruit NA non-pioneer
Rosaceae Prunus myrtifolia (L.) Urb. FALSE FALSE planted diversity non-pioneer
Rubiaceae Genipa americana L. FALSE FALSE planted diversity non-pioneer
Rubiaceae Amaioua intermedia Mart. FALSE TRUE recruit NA non-pioneer
Rubiaceae Palicourea sessilis (Vell.) C.M.Taylor TRUE TRUE recruit NA non-pioneer
Salicaceae Casearia sylvestris Sw. FALSE TRUE recruit NA pioneer
Sapindaceae Matayba elaeagnoides Radlk. FALSE TRUE recruit NA non-pioneer
Sapindaceae Sapindus saponaria L. TRUE TRUE recruit NA non-pioneer
Sapotaceae Pouteria ramiflora Radlk. TRUE TRUE recruit NA non-pioneer
Solanaceae Solanum didymum Dunal FALSE TRUE recruit NA pioneer
Solanaceae Solanum sisymbriifolium Lam. FALSE TRUE recruit NA pioneer
Solanaceae Cestrum corymbosum Schltdl. TRUE FALSE recruit NA non-pioneer
Solanaceae Solanum granuloso-leprosum Dunal TRUE TRUE recruit NA pioneer
Unclassified Unclassified NA TRUE TRUE recruit NA NA
Urticaceae Cecropia pachystachya Trécul TRUE TRUE planted diversity pioneer
Verbenaceae Citharexylum myrianthum Cham. FALSE FALSE planted diversity non-pioneer
Verbenaceae Lantana camara L. TRUE TRUE recruit NA pioneer
Vochysiaceae Qualea grandiflora Mart. FALSE TRUE recruit NA non-pioneer

D) Total AGB

(i) Organizing data

Estimation in Mg.ha-1

576 = 12m x 12m x 4 plots

show code
AGB_data0 <- trees_clean |>
  dplyr::group_by(ano, tratamento, bloco, plantio)|>
  dplyr::summarise(AGB_sum = sum(AGB_ferez, na.rm=T)/1000, .groups = "drop") |> # writexl::write_xlsx("AGB_plot.xlsx")
  dplyr::mutate(
    AGB_Mgha = dplyr::case_when(
      tratamento == "100" ~ (AGB_sum / 576) * 10000,
      tratamento %in% c("F50", "N50") & plantio == "CP" ~ (AGB_sum / 576) * 5000,
      tratamento %in% c("F50", "N50") & plantio == "SP" ~ (AGB_sum / 576) * 5000,
      tratamento %in% c("F25", "N25") & plantio == "CP" ~ (AGB_sum / 576) * 2500,
      tratamento %in% c("F25", "N25") & plantio == "SP" ~ (AGB_sum / 576) * 7500),
    Tratamento = ordered(tratamento, levels = ordem,
                         labels = c("Nuclei 25%", "Strips 25%", "Nuclei 50%", "Strips 50%", "Planting 100%"))) 
AGB_data <- AGB_data0 |> 
  dplyr::group_by(ano, tratamento, bloco)|>
  dplyr::summarise(AGB_Mgha = sum(AGB_Mgha), .groups = "drop") |>
  dplyr::mutate(
    Ano = ano-2021,
    ANO = factor(ano),
    scaled_AGB = scale(AGB_Mgha))

(ii) Analysing Total aboveground biomass per hectare

Comparing models

show code
AGB_model0 <- lme4::lmer(scaled_AGB ~ 0 + (1|bloco), data = AGB_data)
AGB_model1 <- lme4::lmer(scaled_AGB ~ ANO + (1|bloco), data = AGB_data)
AGB_model2 <- lme4::lmer(scaled_AGB ~ ANO + bloco + (1|bloco), data = AGB_data) 
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
unable to evaluate scaled gradient
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Hessian is numerically singular: parameters are not uniquely determined
show code
AGB_model3 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + bloco + (1|bloco), data = AGB_data)   # consider rescaling
AGB_model4 <- lme4::lmer(scaled_AGB ~ tratamento + ANO + bloco + (1|bloco), data = AGB_data) # isSingular
AGB_model5 <- lme4::lmer(scaled_AGB ~ tratamento + ANO + (1|bloco), data = AGB_data)         # best model 
AGB_model6 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + (1|bloco), data = AGB_data) 
AGB_model7 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + (1+Ano|bloco), data = AGB_data) # isSingular
boundary (singular) fit: see help('isSingular')
show code
AIC(AGB_model0, AGB_model1, AGB_model2, AGB_model3, AGB_model4, AGB_model5,AGB_model6, AGB_model7 )
           df       AIC
AGB_model0  2 144.60044
AGB_model1  4 141.90535
AGB_model2  8 144.14482
AGB_model3 16  93.47468
AGB_model4 12  88.21817
AGB_model5  8  85.97871
AGB_model6 12  91.23521
AGB_model7 14  90.49850

Exploring residuals

show code
par(mfrow=c(1,2))
DHARMa::plotResiduals(AGB_model5)
DHARMa::plotQQunif(AGB_model5)

One outlier

Marginal \(R^2\) / Conditional \(R^2\)

show code
MuMIn::r.squaredGLMM(AGB_model5) # without interaction 
           R2m       R2c
[1,] 0.6955092 0.8374889
show code
0.696/ 0.837
[1] 0.8315412
show code
MuMIn::r.squaredGLMM(AGB_model6) # with interaction
          R2m       R2c
[1,] 0.701532 0.8419478
show code
0.702/ 0.842
[1] 0.8337292

0.702/ 0.842 - variance explained by fixed effects and by the entire model (around r0.702/ 0.842 % is explained by fixed effects)

show code
car::Anova(AGB_model6)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: scaled_AGB
                  Chisq Df Pr(>Chisq)    
tratamento     165.4367  4  < 2.2e-16 ***
ANO             47.4491  1  5.645e-12 ***
tratamento:ANO   4.6061  4     0.3302    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

interaction is no significant at all (\(p-value \sim 0.33\))

show code
car::Anova(AGB_model5)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: scaled_AGB
             Chisq Df Pr(>Chisq)    
tratamento 162.968  4  < 2.2e-16 ***
ANO         46.741  1  8.102e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Summary

show code
sjPlot::tab_model(AGB_model5)
  scaled AGB
Predictors Estimates CI p
(Intercept) 0.96 0.51 – 1.42 <0.001
tratamento [F25] -1.83 -2.20 – -1.46 <0.001
tratamento [F50] -1.35 -1.72 – -0.98 <0.001
tratamento [N25] -2.19 -2.56 – -1.82 <0.001
tratamento [N50] -1.44 -1.82 – -1.07 <0.001
ANO [2024] 0.80 0.56 – 1.03 <0.001
Random Effects
σ2 0.17
τ00 bloco 0.15
ICC 0.47
N bloco 5
Observations 50
Marginal R2 / Conditional R2 0.696 / 0.837
show code
broom.mixed::tidy(AGB_model5) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))
# A tibble: 8 × 6
  effect   group    term            estimate std.error statistic
  <chr>    <chr>    <chr>              <dbl>     <dbl>     <dbl>
1 fixed    <NA>     (Intercept)        0.965     0.224      4.31
2 fixed    <NA>     tratamentoF25     -1.83      0.184     -9.92
3 fixed    <NA>     tratamentoF50     -1.35      0.184     -7.34
4 fixed    <NA>     tratamentoN25     -2.19      0.184    -11.9 
5 fixed    <NA>     tratamentoN50     -1.44      0.184     -7.83
6 fixed    <NA>     ANO2024            0.797     0.117      6.84
7 ran_pars bloco    sd__(Intercept)    0.385    NA         NA   
8 ran_pars Residual sd__Observation    0.412    NA         NA   

Assessing comparisons:

show code
emmeans::emmeans(AGB_model6, pairwise ~ tratamento | ANO) |> multcomp::cld()
ANO = 2022:
 tratamento emmean    SE   df lower.CL upper.CL .group
 N25        -1.028 0.251 13.4  -1.5696  -0.4868  1    
 F25        -0.804 0.251 13.4  -1.3450  -0.2622  1    
 N50        -0.527 0.251 13.4  -1.0680   0.0148  1    
 F50        -0.423 0.251 13.4  -0.9645   0.1183  1    
 100         0.790 0.251 13.4   0.2484   1.3312   2   

ANO = 2024:
 tratamento emmean    SE   df lower.CL upper.CL .group
 N25        -0.629 0.251 13.4  -1.1701  -0.0873  1    
 F25        -0.125 0.251 13.4  -0.6663   0.4165  12   
 N50         0.366 0.251 13.4  -0.1756   0.9072   2   
 F50         0.443 0.251 13.4  -0.0981   0.9847   2   
 100         1.936 0.251 13.4   1.3949   2.4777    3  

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 5 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
show code
AGB_model6a <- lme4::lmer(scaled_AGB ~ tratamento*Ano + (1|bloco), data = AGB_data)
AIC(AGB_model6a, AGB_model6)
            df      AIC
AGB_model6a 12 98.16669
AGB_model6  12 91.23521
show code
emmeans::emtrends(AGB_model6a, pairwise ~ tratamento, var = "Ano", poly.trend = TRUE) |> multcomp::cld()
 tratamento Ano.trend    SE df lower.CL upper.CL .group
 N25            0.200 0.129 36  -0.0625    0.462  1    
 F25            0.339 0.129 36   0.0771    0.602  1    
 F50            0.433 0.129 36   0.1709    0.695  1    
 N50            0.446 0.129 36   0.1840    0.708  1    
 100            0.573 0.129 36   0.3110    0.836  1    

Results are averaged over the levels of: Ano 
Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 5 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

(iii) AGB only inside

Preparing data:

Comparison only inside plots over the years

show code
ordem <-c("N25", "F25", "N50", "F50", "100")
AGB_inside0 <- trees_clean |>
  dplyr::group_by(ano, tratamento, bloco, parcela, plantio)|>
  dplyr::summarise(
    AGB_sum = sum(AGB_ferez, na.rm=T),
    .groups = "drop") |>
  dplyr::mutate(
    tratamento = ordered(tratamento, levels = ordem),
    Parcela = paste(tratamento, bloco, parcela, sep="_"),
    Ano = ano-2021) |>
  dplyr::filter(plantio == "CP")

View

show code
AGB_inside0 |>
  ggplot(aes(x = factor(ano), y = log(AGB_sum),  group = Parcela, colour = parcela)) +
  geom_point() +
  geom_line()+facet_grid(bloco~tratamento)+
  theme_minimal()

Plots with decreasing :

show code
deltas <- AGB_inside0 |>
  dplyr::select(ano, Parcela, AGB_sum)|>
  tidyr::pivot_wider(names_from = ano, names_prefix = "ano_", values_from = AGB_sum)|>
  dplyr::mutate(delta_AGB = (ano_2024/ano_2022))|>
  dplyr::arrange(delta_AGB) 
deltas |> head(10)
# A tibble: 10 × 4
   Parcela ano_2022 ano_2024 delta_AGB
   <chr>      <dbl>    <dbl>     <dbl>
 1 100_4_4     422.     216.     0.511
 2 F50_5_3     184.     128.     0.698
 3 N50_5_4     324.     228.     0.701
 4 F25_5_1     269.     198.     0.736
 5 F50_1_1     152.     141.     0.930
 6 N25_3_4     165.     157.     0.954
 7 100_5_3     230.     220.     0.957
 8 F50_3_4     208.     208.     0.998
 9 N50_1_4     485.     493.     1.02 
10 F50_2_4     549.     577.     1.05 
show code
deltas |> tail(10)
# A tibble: 10 × 4
   Parcela ano_2022 ano_2024 delta_AGB
   <chr>      <dbl>    <dbl>     <dbl>
 1 F50_2_1    216.      758.      3.51
 2 N50_5_1    170.      611.      3.58
 3 N50_5_3    153.      554.      3.62
 4 N25_3_2    223.      846.      3.79
 5 F50_2_3     77.5     313.      4.03
 6 F50_2_2     82.8     367.      4.43
 7 F50_3_3    175.      799.      4.56
 8 F25_3_3    114.      557.      4.88
 9 N50_4_2    136.      784.      5.75
10 N50_2_4     88.6     690.      7.79

How much % is planted or recruit in AGB?

show code
trees_clean |>
  dplyr::filter(plantio == "CP") |>
  #dplyr::filter(ano == 2024) |>
  dplyr::group_by(ano, tratamento, bloco)|>
  dplyr::summarise(
    AGB_sum = sum(AGB_ferez, na.rm=T)/1000,
    AGB_planted = sum(AGB_ferez[RP == "planted"], na.rm = TRUE)/1000,
    prop_AGB_planted = (AGB_planted / AGB_sum)*100,
    .groups = "drop") |>  #dplyr::group_by(tratamento,plantio)|>
  dplyr::summarise(mean(prop_AGB_planted),
                   sd(prop_AGB_planted)/sqrt(dplyr::n()))
# A tibble: 1 × 2
  `mean(prop_AGB_planted)` `sd(prop_AGB_planted)/sqrt(dplyr::n())`
                     <dbl>                                   <dbl>
1                     96.9                                   0.758

Should I treat them with a random factor (1|bloco)?

show code
AGB_inside <- trees_clean |>
  dplyr::filter(plantio == "CP") |>
  #dplyr::filter(ano == 2024) |>
  dplyr::group_by(ano, tratamento, bloco)|>
  dplyr::summarise(AGB_sum = sum(AGB_ferez, na.rm=T)/1000, .groups = "drop") |>
  dplyr::mutate(
    Tratamento = ordered(tratamento, levels = ordem),
    Ano = ano-2021,
    ANO = ordered(ano),
    scaled_AGB = scale(AGB_sum)) 

Comparing models

show code
AGB_inside_m0 <- lme4::lmer(scaled_AGB ~ 0 + (1|bloco), data = AGB_inside)
AGB_inside_m1 <- lme4::lmer(scaled_AGB ~ ANO + (1|bloco), data = AGB_inside)
AGB_inside_m2 <- lme4::lmer(scaled_AGB ~ ANO + bloco + (1|bloco), data = AGB_inside)   # best model 
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
unable to evaluate scaled gradient
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Hessian is numerically singular: parameters are not uniquely determined
show code
AGB_inside_m3 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + bloco + (1|bloco), data = AGB_inside)   
AGB_inside_m4 <- lme4::lmer(scaled_AGB ~ tratamento + ANO + bloco + (1|bloco), data = AGB_inside) 
AGB_inside_m5 <- lme4::lmer(scaled_AGB ~ tratamento + ANO + (1|bloco), data = AGB_inside)         
AGB_inside_m6 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + (1|bloco), data = AGB_inside) 
AGB_inside_m7 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + (1+Ano|bloco), data = AGB_inside) # isSingular
boundary (singular) fit: see help('isSingular')
show code
AIC(AGB_inside_m0, AGB_inside_m1, AGB_inside_m2, AGB_inside_m3, AGB_inside_m4,AGB_inside_m5, AGB_inside_m6, AGB_inside_m7)
              df      AIC
AGB_inside_m0  2 138.7732
AGB_inside_m1  4 117.1008
AGB_inside_m2  8 116.0631
AGB_inside_m3 16 122.2545
AGB_inside_m4 12 115.0812
AGB_inside_m5  8 116.1189
AGB_inside_m6 12 123.2922
AGB_inside_m7 14 123.2583

Exploring residuals

show code
par(mfrow=c(1,2))
DHARMa::plotResiduals(AGB_inside_m6)
DHARMa::plotQQunif(AGB_inside_m6)

Marginal \(R^2\) / Conditional \(R^2\)

show code
MuMIn::r.squaredGLMM(AGB_inside_m6)
          R2m      R2c
[1,] 0.396267 0.699652
show code
0.396/  0.6996
[1] 0.5660377

0.396/ 0.6996 variance explained by fixed effects and by the entire model (around 0.5660377 % is explained by fixed effects)

show code
car::Anova(AGB_inside_m6)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: scaled_AGB
                 Chisq Df Pr(>Chisq)    
tratamento     14.3671  4   0.006211 ** 
ANO            47.6142  1  5.189e-12 ***
tratamento:ANO  2.6673  4   0.614946    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

interaction is no significant at all (\(p-value \sim 0.85\))

show code
sjPlot::tab_model(AGB_inside_m5)
  scaled AGB
Predictors Estimates CI p
(Intercept) 0.02 -0.62 – 0.66 0.949
tratamento [F25] 0.50 -0.02 – 1.01 0.059
tratamento [F50] 0.00 -0.51 – 0.52 0.987
tratamento [N25] -0.47 -0.99 – 0.04 0.070
tratamento [N50] -0.13 -0.64 – 0.39 0.622
ANO [linear] 0.80 0.57 – 1.03 <0.001
Random Effects
σ2 0.33
τ00 bloco 0.34
ICC 0.51
N bloco 5
Observations 50
Marginal R2 / Conditional R2 0.390 / 0.702
show code
broom.mixed::tidy(AGB_inside_m5) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))
# A tibble: 8 × 6
  effect   group    term            estimate std.error statistic
  <chr>    <chr>    <chr>              <dbl>     <dbl>     <dbl>
1 fixed    <NA>     (Intercept)        0.02      0.318     0.064
2 fixed    <NA>     tratamentoF25      0.496     0.256     1.94 
3 fixed    <NA>     tratamentoF50      0.004     0.256     0.016
4 fixed    <NA>     tratamentoN25     -0.475     0.256    -1.86 
5 fixed    <NA>     tratamentoN50     -0.127     0.256    -0.497
6 fixed    <NA>     ANO.L              0.802     0.114     7.02 
7 ran_pars bloco    sd__(Intercept)    0.585    NA        NA    
8 ran_pars Residual sd__Observation    0.571    NA        NA    

Assessing comparisons:

show code
emmeans::emmeans(AGB_inside_m6, pairwise ~ tratamento |ANO) |> multcomp::cld()
ANO = 2022:
 tratamento  emmean    SE   df lower.CL upper.CL .group
 N25        -0.8983 0.369 12.2   -1.700  -0.0969  1    
 N50        -0.7449 0.369 12.2   -1.546   0.0564  1    
 F50        -0.5376 0.369 12.2   -1.339   0.2637  1    
 100        -0.3908 0.369 12.2   -1.192   0.4105  1    
 F25        -0.2642 0.369 12.2   -1.066   0.5371  1    

ANO = 2024:
 tratamento  emmean    SE   df lower.CL upper.CL .group
 N25        -0.0107 0.369 12.2   -0.812   0.7906  1    
 100         0.4316 0.369 12.2   -0.370   1.2329  12   
 N50         0.5318 0.369 12.2   -0.269   1.3332  12   
 F50         0.5866 0.369 12.2   -0.215   1.3879  12   
 F25         1.2965 0.369 12.2    0.495   2.0979   2   

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 5 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
show code
AGB_inside_m6a <- lme4::lmer(scaled_AGB ~ tratamento*Ano + (1|bloco), data = AGB_inside) 
AIC(AGB_inside_m6a, AGB_inside_m6)
               df      AIC
AGB_inside_m6a 12 126.7579
AGB_inside_m6  12 123.2922
show code
emmeans::emtrends(AGB_inside_m6a, pairwise ~ tratamento, var = "Ano", poly.trend = TRUE) |> multcomp::cld()
 tratamento Ano.trend    SE df lower.CL upper.CL .group
 100            0.411 0.184 36   0.0385    0.784  1    
 N25            0.444 0.184 36   0.0710    0.817  1    
 F50            0.562 0.184 36   0.1894    0.935  1    
 N50            0.638 0.184 36   0.2656    1.011  1    
 F25            0.780 0.184 36   0.4076    1.153  1    

Results are averaged over the levels of: Ano 
Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 5 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

(iv) AGB only outside (unplanted areas)

Preparing data:

Comparison only outside plots over the years

show code
AGB_outside0 <- trees_clean |>
  dplyr::group_by(ano, tratamento, bloco, parcela, plantio)|>
  dplyr::summarise(
    AGB_sum = sum(AGB_ferez, na.rm=T),
    .groups = "drop") |>
  dplyr::mutate(
    tratamento = ordered(tratamento, levels = ordem),
    Parcela = paste(tratamento, bloco, parcela, sep="_"),
    Ano = ano-2021) |>
  dplyr::filter(plantio == "SP")

View

show code
AGB_outside0 |>
  ggplot(aes(x = factor(ano), y = log(AGB_sum),  group = Parcela, colour = parcela)) +
  geom_point() +
  geom_line()+facet_grid(bloco~tratamento)+
  theme_minimal()

it is a mess

Tratar os valores como (1|bloco)?

show code
AGB_outside <- trees_clean |>
  dplyr::filter(plantio == "SP") |>
  #dplyr::filter(ano == 2024) |>
  dplyr::group_by(ano, tratamento, bloco)|>
  dplyr::summarise(AGB_sum = sum(AGB_ferez, na.rm=T)/1000, .groups = "drop") |>
  dplyr::mutate(
    Tratamento = ordered(tratamento, levels = ordem),
    Ano = ano-2021,
    ANO = ordered(ano),
    scaled_AGB = scale(AGB_sum)) 

Comparing models

show code
AGB_outside_m0 <- lme4::lmer(scaled_AGB ~ 0 + (1|bloco), data = AGB_outside) # best model
boundary (singular) fit: see help('isSingular')
show code
AGB_outside_m1 <- lme4::lmer(scaled_AGB ~ ANO + (1|bloco), data = AGB_outside)
AGB_outside_m2 <- lme4::lmer(scaled_AGB ~ ANO + bloco + (1|bloco), data = AGB_outside)   # best model 
AGB_outside_m3 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + bloco + (1|bloco), data = AGB_outside)   
AGB_outside_m4 <- lme4::lmer(scaled_AGB ~ tratamento + ANO + bloco + (1|bloco), data = AGB_outside) 
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
unable to evaluate scaled gradient
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Hessian is numerically singular: parameters are not uniquely determined
show code
AGB_outside_m5 <- lme4::lmer(scaled_AGB ~ tratamento + ANO + (1|bloco), data = AGB_outside)         
AGB_outside_m6 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + (1|bloco), data = AGB_outside) 
AGB_outside_m7 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + (1+Ano|bloco), data = AGB_outside) # isSingular
boundary (singular) fit: see help('isSingular')
show code
AIC(AGB_outside_m0, AGB_outside_m1, AGB_outside_m2, AGB_outside_m3, 
    AGB_outside_m4,AGB_outside_m5, AGB_outside_m6, AGB_outside_m7)
               df      AIC
AGB_outside_m0  2 116.5024
AGB_outside_m1  4 117.0739
AGB_outside_m2  8 120.3459
AGB_outside_m3 14 128.0330
AGB_outside_m4 11 126.0844
AGB_outside_m5  7 122.8125
AGB_outside_m6 10 124.7610
AGB_outside_m7 12 126.6309

Exploring residuals

show code
par(mfrow=c(1,2))
DHARMa::plotResiduals(AGB_outside_m1)
DHARMa::plotQQunif(AGB_outside_m0)

Marginal \(R^2\) / Conditional \(R^2\)

show code
MuMIn::r.squaredGLMM(AGB_outside_m1)
           R2m       R2c
[1,] 0.1505529 0.1714753
show code
0.1505529 /0.1714753
[1] 0.8779859

0.1505529 /0.1714753 variance explained by fixed effects and by the entire model (around 0.8779859 % is explained by fixed effects)

show code
car::Anova(AGB_outside_m5)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: scaled_AGB
            Chisq Df Pr(>Chisq)  
tratamento 0.6275  3    0.89012  
ANO        6.5922  1    0.01024 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

interaction is no significant at all (\(p-value \sim 0.85\))

show code
sjPlot::tab_model(AGB_outside_m5)
  scaled AGB
Predictors Estimates CI p
(Intercept) 0.20 -0.42 – 0.83 0.514
tratamento [F50] -0.30 -1.17 – 0.57 0.493
tratamento [N25] -0.23 -1.10 – 0.64 0.597
tratamento [N50] -0.29 -1.15 – 0.58 0.509
ANO [linear] 0.55 0.11 – 0.98 0.015
Random Effects
σ2 0.91
τ00 bloco 0.01
ICC 0.01
N bloco 5
Observations 40
Marginal R2 / Conditional R2 0.154 / 0.167
show code
broom.mixed::tidy(AGB_outside_m5) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))
# A tibble: 7 × 6
  effect   group    term            estimate std.error statistic
  <chr>    <chr>    <chr>              <dbl>     <dbl>     <dbl>
1 fixed    <NA>     (Intercept)        0.202     0.307     0.66 
2 fixed    <NA>     tratamentoF50     -0.296     0.427    -0.694
3 fixed    <NA>     tratamentoN25     -0.228     0.427    -0.534
4 fixed    <NA>     tratamentoN50     -0.285     0.427    -0.667
5 fixed    <NA>     ANO.L              0.549     0.214     2.57 
6 ran_pars bloco    sd__(Intercept)    0.116    NA        NA    
7 ran_pars Residual sd__Observation    0.956    NA        NA    

Assessing comparisons:

show code
emmeans::emmeans(AGB_outside_m6, pairwise ~ tratamento |ANO) |> multcomp::cld()
ANO = 2022:
 tratamento  emmean    SE df lower.CL upper.CL .group
 F50        -0.6400 0.436 32  -1.5273    0.247  1    
 N25        -0.4221 0.436 32  -1.3094    0.465  1    
 F25        -0.3904 0.436 32  -1.2776    0.497  1    
 N50        -0.0995 0.436 32  -0.9868    0.788  1    

ANO = 2024:
 tratamento  emmean    SE df lower.CL upper.CL .group
 N50        -0.0659 0.436 32  -0.9531    0.821  1    
 N25         0.3708 0.436 32  -0.5165    1.258  1    
 F50         0.4519 0.436 32  -0.4354    1.339  1    
 F25         0.7952 0.436 32  -0.0921    1.682  1    

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 4 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Table S2

Only inside

show code
AGB_inside |> 
  dplyr::group_by(ano, tratamento) |>
  dplyr::summarise(
    Mean = round(mean(AGB_sum),2),
    SE = round(sd(AGB_sum)/sqrt(dplyr::n()), 2),
    .groups="drop") |>
  dplyr::mutate(values = paste(Mean, "pm", SE)) |>
  dplyr::select(-c(Mean, SE)) |>
  tidyr::pivot_wider(names_from = tratamento, values_from = values)
# A tibble: 2 × 6
    ano `100`        F25         F50          N25          N50         
  <dbl> <chr>        <chr>       <chr>        <chr>        <chr>       
1  2022 1.23 pm 0.15 1.32 pm 0.2 1.12 pm 0.26 0.85 pm 0.2  0.96 pm 0.12
2  2024 1.85 pm 0.18 2.5 pm 0.37 1.96 pm 0.46 1.52 pm 0.37 1.92 pm 0.24

Increment per year:

show code
AGB_inside |>
  dplyr::group_by(tratamento) |>
  dplyr::summarise(
    variacao_media = mean((AGB_sum[ano == 2024] - AGB_sum[ano == 2022])/2),
    desvio_padrao = sd((AGB_sum[ano == 2024] - AGB_sum[ano == 2022])/2),
    erro_padrao = sd((AGB_sum[ano == 2024] - AGB_sum[ano == 2022])/2) / sqrt(dplyr::n()),
    min_variacao = min((AGB_sum[ano == 2024] - AGB_sum[ano == 2022])/2),
    max_variacao = max((AGB_sum[ano == 2024] - AGB_sum[ano == 2022])/2),
    .groups="drop")|>
  dplyr::mutate(slope = paste(round(variacao_media,2), "pm", round(erro_padrao,2)))|>
  dplyr::select(tratamento, slope)
# A tibble: 5 × 2
  tratamento slope       
  <chr>      <chr>       
1 100        0.31 pm 0.03
2 F25        0.59 pm 0.08
3 F50        0.42 pm 0.08
4 N25        0.33 pm 0.06
5 N50        0.48 pm 0.1 

Total Aboveground biomass

show code
tabS1_agb_total <- AGB_data |> 
  dplyr::group_by(ano, tratamento)|>
  dplyr::summarise(Mean = round(mean(AGB_Mgha),2),
                   SE = round(sd(AGB_Mgha)/sqrt(dplyr::n()), 2),
                   .groups="drop") 
tabS1_agb_total |>
  dplyr::mutate(values = paste(Mean, "pm", SE)) |>
  dplyr::select(-c(Mean, SE)) |>
  tidyr::pivot_wider(names_from = tratamento, values_from = values)
# A tibble: 2 × 6
    ano `100`         F25           F50           N25          N50          
  <dbl> <chr>         <chr>         <chr>         <chr>        <chr>        
1  2022 21.35 pm 2.64 6.47 pm 0.92  10.03 pm 2.29 4.37 pm 1.07 9.06 pm 0.95 
2  2024 32.06 pm 3.11 12.81 pm 2.12 18.12 pm 4.4  8.11 pm 1.64 17.39 pm 1.95

Increment per year:

show code
AGB_data |>
  dplyr::group_by(tratamento) %>%
  dplyr::summarise(
    variacao_media = mean((AGB_Mgha[ano == 2024] - AGB_Mgha[ano == 2022])/2),
    desvio_padrao = sd((AGB_Mgha[ano == 2024] - AGB_Mgha[ano == 2022])/2),
    erro_padrao = sd((AGB_Mgha[ano == 2024] - AGB_Mgha[ano == 2022])/2) / sqrt(dplyr::n()),
    min_variacao = min((AGB_Mgha[ano == 2024] - AGB_Mgha[ano == 2022])/2),
    max_variacao = max((AGB_Mgha[ano == 2024] - AGB_Mgha[ano == 2022])/2),
    .groups="drop")|>
  dplyr::mutate(slope = paste(round(variacao_media,2), "pm", round(erro_padrao,2)))|>
  dplyr::select(tratamento, slope)
# A tibble: 5 × 2
  tratamento slope       
  <chr>      <chr>       
1 100        5.35 pm 0.59
2 F25        3.17 pm 0.51
3 F50        4.05 pm 0.82
4 N25        1.87 pm 0.35
5 N50        4.17 pm 0.9 

Figure Biomass

Percentuals

show code
agb_plot_data <-AGB_data0 |> 
  dplyr::group_by(ano, Tratamento, plantio, bloco)|>
  dplyr::summarise(
    AGB_Mgha = sum(AGB_Mgha),
    .groups = "drop") |>
 dplyr::group_by(ano, Tratamento, plantio)|>
    dplyr::summarise(meanAGB = mean(AGB_Mgha),
    sd = sd(AGB_Mgha),
    n = dplyr::n(),
    se = sd/sqrt(n),
    .groups = "drop") |>
  dplyr::mutate(
    Ano = as.factor(ifelse(ano == 2022, "5", "7")),
    Tratamento = forcats::fct_recode(Tratamento,
                                     "Nuclei\n25%" = "Nuclei 25%",
                                     "Strips\n25%" = "Strips 25%",
                                     "Nuclei\n50%" = "Nuclei 50%",
                                     "Strips\n50%" = "Strips 50%",
                                     "Planting\n100%" = "Planting 100%"
                                     ))
max_year <- agb_plot_data |>
  dplyr::group_by(ano, Tratamento)|>
  dplyr::summarise(meanAGB = max(meanAGB),.groups="drop")|>
  dplyr::group_by(ano)|>
  dplyr::summarise(max_year = max(meanAGB), .groups="drop")
agb_plot_data |>
  dplyr::group_by(ano, Tratamento)|>
  dplyr::summarise(meanAGB = max(meanAGB),.groups="drop")|>
  dplyr::right_join(max_year, by="ano") |>
  dplyr::mutate(porcentos = round( (meanAGB/max_year)*100, 0))
# A tibble: 2 × 5
    ano Tratamento meanAGB max_year porcentos
  <dbl> <ord>        <dbl>    <dbl>     <dbl>
1  2022 <NA>          48.9     48.9       100
2  2024 <NA>          83.2     83.2       100

Plotting

show code
error_bars <- agb_plot_data |> 
  arrange(Ano, Tratamento, desc(plantio)) |> 
  group_by(Ano, Tratamento) |>
  mutate(meanAGB_new = cumsum(meanAGB)) |>
  ungroup()
porcentos <- data.frame(
  Ano = rep(c("5", "7"), each = 5),  
  Tratamento = rep(unique(agb_plot_data$Tratamento), times = length(unique(agb_plot_data$Ano))),
  meanAGB = c(8, 9, 12, 15, 27, 12, 17, 23, 25, 38),
  label = c(c("17%","27%", "39%", "46%", "100%",
              "21%", "34%", "52%", "53%", "100%"))
  )
text_data <- data.frame(
  Ano = rep(c("5", "7"), each = 5),  
  Tratamento = rep(unique(agb_plot_data$Tratamento), times = length(unique(agb_plot_data$Ano))),
  meanAGB = c(10, 11, 14, 17, 29, 14, 19, 25, 27, 40),
  label = c(c("a","a", "a", "a", "b",
              "A", "AB", "B", "B", "C"))
  )
fig_biomass<-agb_plot_data |>
  ggplot(aes(x = Ano, y = meanAGB, fill = plantio)) +
  geom_bar(stat = "identity", color = "black") +
  geom_errorbar(data=error_bars,
                aes(ymin = meanAGB_new - se, ymax = meanAGB_new + se),
                width = .3) +
  scale_fill_manual(
    labels = c('Planted', 'Unplanted'), 
    values = c("#92D050", "white")
    ) +
  labs(fill = "", 
       y =  expression("Total aboveground biomass (Mg.ha"^{-1} ~")"),
       x = "") +
  facet_grid(.~Tratamento, switch = "both")+
  geom_text(data = text_data, aes(x = Ano, y = meanAGB, label = label), size = 3, color = "black", inherit.aes = FALSE) +
  geom_text(data = porcentos, aes(x = Ano, y = meanAGB, label = label), size = 3, color = "black", inherit.aes = FALSE) +
  scale_y_continuous(expand = c(0, 0), limits = c(0,43))+
  ggthemes::theme_clean()+
  theme(
    #plot.margin = margin(b = 0, l = 0),
    panel.grid.major.y = element_blank(),
    panel.grid.major.x = element_blank(),
    axis.text.x = element_text(color="black", hjust = 1, size=10),
    axis.text.y = element_text(color="black", size=9),
    axis.title.y = element_text(color="black", size=11),
    strip.text = element_text(color="black", hjust = 0.5, size=10),
    panel.spacing = unit(0.1, "mm"),
    strip.placement = "outside",
    legend.text = element_text(size = 11),
    legend.title = element_blank(),
    legend.key.size = unit(3, "mm"),
    legend.box.background = element_blank(),
    legend.box.margin = margin(t = 5, r = 0, b = 0, l = 0, "mm"),
    legend.background = element_blank(), 
    legend.margin = margin(0,0,0,0),
    legend.position = c(0.18, 0.90)
  )
fig_biomass
Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_bar()`).

show code
# ggsave("fig_biomass.png", plot=fig_biomass, width=90, height = 90 ,unit="mm")

II. Results - LiDAR

Canopy Height Metric (CHM_mean, CHM_SD, and CHM_CV), Canopy cover (Canopy_cover), Rugosity (Rugosity_SD, and Rugosity_CV), Leaf Area Index (LAI_mean, LAI_SD, and LAI_CV), and LAIunder (LAIunder_mean, LAIunder_SD, LAIunder_CV).

a) Canopy Height Model

2022

show code
CHM2022 <- lidar_wide |>
  dplyr::filter(ano==2022)|>
  with(lme4::lmer(CHM_mean ~ Tratamento + (1|Bloco) ) )
par(mfrow=c(1,2))
DHARMa::plotResiduals(CHM2022)
DHARMa::plotQQunif(CHM2022)

2024

show code
CHM2024 <- lidar_wide |>
  dplyr::filter(ano==2024)|>
  with(lme4::lmer(CHM_mean ~ Tratamento + (1|Bloco) ) )
par(mfrow=c(1,2))
DHARMa::plotResiduals(CHM2024)
DHARMa::plotQQunif(CHM2024)

Together

show code
modeloCHM <- lidar_wide |>
  with(lme4::lmer(CHM_mean ~ Tratamento*ano + (1|Bloco) ) )
par(mfrow=c(1,2))
DHARMa::plotResiduals(modeloCHM)
DHARMa::plotQQunif(modeloCHM)

Comparing estimated means

show code
emmeans::emmeans(modeloCHM, pairwise ~ Tratamento | ano) |> multcomp::cld()
ano = 2022:
 Tratamento    emmean    SE df lower.CL upper.CL .group
 Nuclei 25%      1.89 0.437 10    0.919     2.87  1    
 Strips 25%      2.12 0.437 10    1.142     3.09  1    
 Strips 50%      2.93 0.437 10    1.956     3.90  1    
 Nuclei 50%      3.00 0.437 10    2.027     3.97  1    
 Planting 100%   4.91 0.437 10    3.941     5.89   2   

ano = 2024:
 Tratamento    emmean    SE df lower.CL upper.CL .group
 Nuclei 25%      2.76 0.437 10    1.788     3.73  1    
 Strips 25%      3.16 0.437 10    2.187     4.13  1    
 Strips 50%      4.43 0.437 10    3.454     5.40   2   
 Nuclei 50%      4.70 0.437 10    3.724     5.67   2   
 Planting 100%   6.99 0.437 10    6.018     7.96    3  

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 5 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Comparing estimated slopes

show code
emmeans::emtrends(modeloCHM, pairwise ~ Tratamento, var = "ano", poly.trend = TRUE)|> multcomp::cld()
 Tratamento    ano.trend    SE df lower.CL upper.CL .group
 Nuclei 25%        0.434 0.201 36    0.026    0.843  1    
 Strips 25%        0.523 0.201 36    0.114    0.931  1    
 Strips 50%        0.749 0.201 36    0.341    1.157  1    
 Nuclei 50%        0.849 0.201 36    0.440    1.257  1    
 Planting 100%     1.039 0.201 36    0.630    1.447  1    

Results are averaged over the levels of: ano 
Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 5 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Comparing fixed and random factors variance

show code
MuMIn::r.squaredGLMM(modeloCHM)
           R2m       R2c
[1,] 0.7002626 0.8728069
show code
0.7002626  /0.8728069
[1] 0.802311

Results:

The trend of variation around the mean CHM was the same over time among treatments (n.s.). The significant differences (\(p<0.001\)) observed between the 100% planting treatment and the other treatments in 2022 persisted in 2024.

show code
sjPlot::tab_model(modeloCHM)
  Dependent variable
Predictors Estimates CI p
(Intercept) -876.29 -1700.80 – -51.79 0.038
Tratamento [Strips 25%] -178.32 -1344.35 – 987.70 0.759
Tratamento [Nuclei 50%] -836.96 -2002.98 – 329.07 0.154
Tratamento [Strips 50%] -635.54 -1801.57 – 530.48 0.277
Tratamento [Planting
100%]
-1218.82 -2384.85 – -52.80 0.041
ano 0.43 0.03 – 0.84 0.037
Tratamento [Strips 25%] ×
ano
0.09 -0.49 – 0.66 0.758
Tratamento [Nuclei 50%] ×
ano
0.41 -0.16 – 0.99 0.154
Tratamento [Strips 50%] ×
ano
0.31 -0.26 – 0.89 0.276
Tratamento [Planting
100%] × ano
0.60 0.03 – 1.18 0.040
Random Effects
σ2 0.41
τ00 Bloco 0.55
ICC 0.58
N Bloco 5
Observations 50
Marginal R2 / Conditional R2 0.700 / 0.873
show code
broom.mixed::tidy(modeloCHM) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))
# A tibble: 12 × 6
   effect   group    term                         estimate std.error statistic
   <chr>    <chr>    <chr>                           <dbl>     <dbl>     <dbl>
 1 fixed    <NA>     (Intercept)                  -876.      407.        -2.15
 2 fixed    <NA>     TratamentoStrips 25%         -178.      576.        -0.31
 3 fixed    <NA>     TratamentoNuclei 50%         -837.      576.        -1.45
 4 fixed    <NA>     TratamentoStrips 50%         -636.      576.        -1.10
 5 fixed    <NA>     TratamentoPlanting 100%     -1219.      576.        -2.12
 6 fixed    <NA>     ano                             0.434     0.201      2.16
 7 fixed    <NA>     TratamentoStrips 25%:ano        0.088     0.285      0.31
 8 fixed    <NA>     TratamentoNuclei 50%:ano        0.414     0.285      1.46
 9 fixed    <NA>     TratamentoStrips 50%:ano        0.315     0.285      1.11
10 fixed    <NA>     TratamentoPlanting 100%:ano     0.604     0.285      2.12
11 ran_pars Bloco    sd__(Intercept)                 0.742    NA         NA   
12 ran_pars Residual sd__Observation                 0.637    NA         NA   

Table S2

Values per year

show code
lidar_wide |>
  dplyr::group_by(ano, Tratamento)|>
  dplyr::summarise(Mean = round(mean(CHM_mean), 2),
                   SE = round(sd(CHM_mean)/sqrt(dplyr::n()), 2)) |>
  dplyr::mutate(values = paste(Mean, "pm", SE)) |>
  dplyr::select(-c(Mean, SE)) |>
  tidyr::pivot_wider(names_from = Tratamento, values_from = values)
`summarise()` has grouped output by 'ano'. You can override using the `.groups`
argument.
# A tibble: 2 × 6
# Groups:   ano [2]
    ano `Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%` `Planting 100%`
  <dbl> <chr>        <chr>        <chr>        <chr>        <chr>          
1  2022 1.89 pm 0.3  2.12 pm 0.41 3 pm 0.26    2.93 pm 0.46 4.91 pm 0.33   
2  2024 2.76 pm 0.48 3.16 pm 0.56 4.7 pm 0.28  4.43 pm 0.69 6.99 pm 0.4    

“Slope” (increment/year)

show code
lidar_wide |>
  dplyr::group_by(Tratamento, ano) |>
  dplyr::summarise(
    media = mean(CHM_mean, na.rm = TRUE),
    se = sd(CHM_mean, na.rm = TRUE) / sqrt(dplyr::n()),
    .groups = "drop"
  ) |>
  tidyr::pivot_wider(
    names_from = ano,
    values_from = c(media, se),
    names_sep = "_"
  ) |>
  dplyr::mutate(
    variacao_media = round((media_2024 - media_2022) / 2, 2),
    SE_variacao = round(sqrt(se_2022^2 + se_2024^2) / 2, 2)) |>
  dplyr::mutate(values = paste(variacao_media, "pm", SE_variacao)) |>
  dplyr::select(Tratamento, values) |>
  tidyr::pivot_wider(names_from = Tratamento, values_from = values)
# A tibble: 1 × 5
  `Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%` `Planting 100%`
  <chr>        <chr>        <chr>        <chr>        <chr>          
1 0.43 pm 0.28 0.52 pm 0.35 0.85 pm 0.19 0.75 pm 0.42 1.04 pm 0.26   

b) Canopy cover

show code
lidar_wide |>
  dplyr::group_by(ano)|>
  dplyr::summarise(mean(Canopy_cover),
                   min(Canopy_cover),
                   max(Canopy_cover),
                   sd(Canopy_cover))
# A tibble: 2 × 5
    ano `mean(Canopy_cover)` `min(Canopy_cover)` `max(Canopy_cover)`
  <dbl>                <dbl>               <dbl>               <dbl>
1  2022                 38.4                5.29                99.8
2  2024                 67.6               10.7                 99.8
# ℹ 1 more variable: `sd(Canopy_cover)` <dbl>

2022

show code
CC2022 <- lidar_wide |>
  dplyr::filter(ano==2022)|>
  with( glmmTMB::glmmTMB(Canopy_cover ~ Tratamento + (1|Bloco)) )
Warning in glmmTMB::glmmTMB(Canopy_cover ~ Tratamento + (1 | Bloco)): use of
the 'data' argument is recommended
show code
CC2022a <- lidar_wide |>
  dplyr::filter(ano==2022)|>
  with(glmmTMB::glmmTMB(Canopy_cover/100 ~ Tratamento + (1 | Bloco),
                        family = glmmTMB::beta_family()))
Warning in glmmTMB::glmmTMB(Canopy_cover/100 ~ Tratamento + (1 | Bloco), : use
of the 'data' argument is recommended
show code
AIC(CC2022, CC2022a)
        df       AIC
CC2022   7 228.16951
CC2022a  7  -5.60719
show code
par(mfrow=c(2,2))
DHARMa::plotResiduals(CC2022a)
DHARMa::plotQQunif(CC2022a)
DHARMa::plotResiduals(CC2022)
DHARMa::plotQQunif(CC2022)

2024

show code
CC2024 <- lidar_wide |>
  dplyr::filter(ano==2024) |>
  with(lme4::lmer(Canopy_cover ~ Tratamento + (1|Bloco) ) )
CC2024a <- lidar_wide |>
  dplyr::filter(ano==2024) |>
  with(glmmTMB::glmmTMB(Canopy_cover/100 ~ Tratamento + (1 | Bloco),
                        family = glmmTMB::beta_family()))
Warning in glmmTMB::glmmTMB(Canopy_cover/100 ~ Tratamento + (1 | Bloco), : use
of the 'data' argument is recommended
show code
AIC(CC2024, CC2024a)
        df       AIC
CC2024   7 212.04388
CC2024a  7 -16.11127
show code
par(mfrow=c(2,2))
DHARMa::plotResiduals(CC2024a)
DHARMa::plotQQunif(CC2024a)
DHARMa::plotResiduals(CC2024)
DHARMa::plotQQunif(CC2024)

models

show code
modeloCC <- lidar_wide |>
  with(lme4::lmer(Canopy_cover ~ Tratamento*ano + (1|Bloco) ) )

modeloCCa <- lidar_wide |>
  with(glmmTMB::glmmTMB(Canopy_cover/100 ~ Tratamento*ano + (1 | Bloco),
                        family = glmmTMB::beta_family()))
Warning in glmmTMB::glmmTMB(Canopy_cover/100 ~ Tratamento * ano + (1 | Bloco),
: use of the 'data' argument is recommended
Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
NA/NaN function evaluation
Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
NA/NaN function evaluation
show code
par(mfrow=c(2,2))
DHARMa::plotResiduals(modeloCC)
DHARMa::plotQQunif(modeloCC)
DHARMa::plotResiduals(modeloCCa)
DHARMa::plotQQunif(modeloCCa)

1 outlier

Estimated means

show code
emmeans::emmeans(modeloCC, pairwise ~ Tratamento | ano) |> multcomp::cld()
ano = 2022:
 Tratamento    emmean   SE   df lower.CL upper.CL .group
 Strips 25%      23.7 11.5 20.8    -0.25     47.6  1    
 Nuclei 50%      26.5 11.5 20.8     2.56     50.4  1    
 Strips 50%      31.3 11.5 20.8     7.42     55.3  1    
 Nuclei 25%      31.6 11.5 20.8     7.70     55.5  1    
 Planting 100%   79.1 11.5 20.8    55.14    103.0   2   

ano = 2024:
 Tratamento    emmean   SE   df lower.CL upper.CL .group
 Strips 25%      44.7 11.5 20.8    20.78     68.6  1    
 Nuclei 25%      56.5 11.5 20.8    32.62     80.5  1    
 Strips 50%      69.8 11.5 20.8    45.88     93.7  12   
 Nuclei 50%      70.9 11.5 20.8    47.02     94.9  12   
 Planting 100%   96.2 11.5 20.8    72.29    120.1   2   

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 5 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Estimated slopes

show code
emmeans::emtrends(modeloCC, pairwise ~ Tratamento, var = "ano", poly.trend = TRUE)|> multcomp::cld()
 Tratamento    ano.trend   SE df lower.CL upper.CL .group
 Planting 100%      8.57 6.71 36    -5.02     22.2  1    
 Strips 25%        10.52 6.71 36    -3.08     24.1  1    
 Nuclei 25%        12.46 6.71 36    -1.14     26.1  1    
 Strips 50%        19.23 6.71 36     5.63     32.8  1    
 Nuclei 50%        22.23 6.71 36     8.63     35.8  1    

Results are averaged over the levels of: ano 
Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 5 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Comparing random and fixed effects

show code
MuMIn::r.squaredGLMM(modeloCC)
          R2m       R2c
[1,] 0.469023 0.6387404
show code
0.469023 / 0.6387404
[1] 0.7342936
show code
sjPlot::tab_model(modeloCC)
  Dependent variable
Predictors Estimates CI p
(Intercept) -25156.96 -52617.32 – 2303.39 0.071
Tratamento [Strips 25%] 3917.03 -34917.77 – 42751.83 0.839
Tratamento [Nuclei 50%] -19762.19 -58596.99 – 19072.61 0.309
Tratamento [Strips 50%] -13690.08 -52524.88 – 25144.72 0.480
Tratamento [Planting
100%]
7898.75 -30936.05 – 46733.55 0.683
ano 12.46 -1.12 – 26.03 0.071
Tratamento [Strips 25%] ×
ano
-1.94 -21.14 – 17.26 0.839
Tratamento [Nuclei 50%] ×
ano
9.77 -9.43 – 28.97 0.309
Tratamento [Strips 50%] ×
ano
6.77 -12.43 – 25.97 0.480
Tratamento [Planting
100%] × ano
-3.88 -23.08 – 15.31 0.684
Random Effects
σ2 449.60
τ00 Bloco 211.22
ICC 0.32
N Bloco 5
Observations 50
Marginal R2 / Conditional R2 0.469 / 0.639
show code
broom.mixed::tidy(modeloCC) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))
# A tibble: 12 × 6
   effect   group    term                         estimate std.error statistic
   <chr>    <chr>    <chr>                           <dbl>     <dbl>     <dbl>
 1 fixed    <NA>     (Intercept)                 -25157.    13565.      -1.86 
 2 fixed    <NA>     TratamentoStrips 25%          3917.    19183.       0.204
 3 fixed    <NA>     TratamentoNuclei 50%        -19762.    19183.      -1.03 
 4 fixed    <NA>     TratamentoStrips 50%        -13690.    19183.      -0.714
 5 fixed    <NA>     TratamentoPlanting 100%       7899.    19183.       0.412
 6 fixed    <NA>     ano                             12.5       6.70     1.86 
 7 fixed    <NA>     TratamentoStrips 25%:ano        -1.94      9.48    -0.205
 8 fixed    <NA>     TratamentoNuclei 50%:ano         9.77      9.48     1.03 
 9 fixed    <NA>     TratamentoStrips 50%:ano         6.77      9.48     0.714
10 fixed    <NA>     TratamentoPlanting 100%:ano     -3.88      9.48    -0.409
11 ran_pars Bloco    sd__(Intercept)                 14.5      NA       NA    
12 ran_pars Residual sd__Observation                 21.2      NA       NA    

Table S2

per year

show code
lidar_wide |>
  dplyr::group_by(ano, Tratamento)|>
  dplyr::summarise(Mean = round(mean(Canopy_cover), 2),
                   SE = round((sd(Canopy_cover)/sqrt(dplyr::n())), 2),
                   .groups = "drop") |>
  dplyr::mutate(values = paste(Mean, "pm", SE)) |>
  dplyr::select(-c(Mean, SE)) |>
  tidyr::pivot_wider(names_from = Tratamento, values_from = values)
# A tibble: 2 × 6
    ano `Nuclei 25%`   `Strips 25%`  `Nuclei 50%`  `Strips 50%`  `Planting 100%`
  <dbl> <chr>          <chr>         <chr>         <chr>         <chr>          
1  2022 31.62 pm 17.23 23.67 pm 4.99 26.48 pm 1.81 31.34 pm 7.74 79.06 pm 10.54 
2  2024 56.54 pm 18.23 44.7 pm 13.14 70.94 pm 11.3 69.8 pm 13.76 96.2 pm 1.89   

“Slope” (increment/year)

show code
lidar_wide |>
  dplyr::group_by(Tratamento, ano) |>
  dplyr::summarise(
    media = mean(Canopy_cover, na.rm = TRUE),
    se = sd(Canopy_cover, na.rm = TRUE) / sqrt(dplyr::n()),
    .groups = "drop"
  ) |>
  tidyr::pivot_wider(
    names_from = ano,
    values_from = c(media, se),
    names_sep = "_"
  ) |>
  dplyr::mutate(
    variacao_media = round((media_2024 - media_2022) / 2, 2),
    SE_variacao = round(sqrt(se_2022^2 + se_2024^2) / 2, 2)) |>
  dplyr::mutate(values = paste(variacao_media, "pm", SE_variacao)) |>
  dplyr::select(Tratamento, values) |>
  tidyr::pivot_wider(names_from = Tratamento, values_from = values)
# A tibble: 1 × 5
  `Nuclei 25%`   `Strips 25%`  `Nuclei 50%`  `Strips 50%`  `Planting 100%`
  <chr>          <chr>         <chr>         <chr>         <chr>          
1 12.46 pm 12.54 10.52 pm 7.03 22.23 pm 5.72 19.23 pm 7.89 8.57 pm 5.35   

c) variance of Canopy cover

Performing a pairwise comparison of variances is not as straightforward as comparing means, but you can use specific statistical techniques to compare variances between treatment groups. Here’s how you can do this in R:

show code
pairwise_levene <- function(data, response, group) {
  treatment_pairs <- combn(unique(data[[group]]), 2, simplify = FALSE)
  p_values <- vector("list", length(treatment_pairs))
  
  for (i in seq_along(treatment_pairs)) { #i =1
    subset_data <- data[data[[group]] %in% treatment_pairs[[i]], ]
    test_result <- leveneTest(subset_data[[response]] ~ subset_data[[group]])
    p_values[[i]] <- test_result[1, "Pr(>F)"]
  }
  
  comparison <- sapply(treatment_pairs, function(x) paste(x, collapse = " vs "))
  results <- data.frame(Comparison = comparison, P_Value = unlist(p_values))
  results$Adjusted_P_Value <- p.adjust(results$P_Value, method = "bonferroni")
  
  return(results)
}

pairwise_levene(lidar2022_wide , "Canopy_cover", "Tratamento")
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
                     Comparison    P_Value Adjusted_P_Value
1      25_faixa vs 25_tabuleiro 0.41911844        1.0000000
2          25_faixa vs 50_faixa 0.60465350        1.0000000
3      25_faixa vs 50_tabuleiro 0.52704716        1.0000000
4      25_faixa vs 100_controle 0.21179066        1.0000000
5      25_tabuleiro vs 50_faixa 0.56742080        1.0000000
6  25_tabuleiro vs 50_tabuleiro 0.31075886        1.0000000
7  25_tabuleiro vs 100_controle 0.89320378        1.0000000
8      50_faixa vs 50_tabuleiro 0.26004211        1.0000000
9      50_faixa vs 100_controle 0.43307531        1.0000000
10 50_tabuleiro vs 100_controle 0.07724271        0.7724271
show code
pairwise_levene(lidar2024_wide , "Canopy_cover", "Tratamento")
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
                     Comparison    P_Value Adjusted_P_Value
1      25_faixa vs 25_tabuleiro 0.47523130        1.0000000
2          25_faixa vs 50_faixa 0.73908327        1.0000000
3      25_faixa vs 50_tabuleiro 0.87724043        1.0000000
4      25_faixa vs 100_controle 0.28751774        1.0000000
5      25_tabuleiro vs 50_faixa 0.65306786        1.0000000
6  25_tabuleiro vs 50_tabuleiro 0.49847386        1.0000000
7  25_tabuleiro vs 100_controle 0.08104640        0.8104640
8      50_faixa vs 50_tabuleiro 0.81722803        1.0000000
9      50_faixa vs 100_controle 0.11052807        1.0000000
10 50_tabuleiro vs 100_controle 0.09110862        0.9110862

Table S2 (Rugosity)

show code
lidar_wide |>
  dplyr::group_by(ano, Tratamento)|>
  dplyr::summarise(Mean = round(mean(Rugosity_SD),2),
                   SE = round((sd(Rugosity_SD)/sqrt(dplyr::n())), 2),
                   .groups="drop") |>
  dplyr::mutate(values = paste(Mean, "pm", SE)) |>
  dplyr::select(-c(Mean, SE)) |>
  tidyr::pivot_wider(names_from = Tratamento, values_from = values)
# A tibble: 2 × 6
    ano `Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%` `Planting 100%`
  <dbl> <chr>        <chr>        <chr>        <chr>        <chr>          
1  2022 2.08 pm 0.22 2.38 pm 0.17 2.73 pm 0.31 2.64 pm 0.19 2.57 pm 0.06   
2  2024 2.77 pm 0.32 3.15 pm 0.18 3.41 pm 0.36 3.3 pm 0.17  2.79 pm 0.09   

Increment/year

show code
lidar_wide |>
  dplyr::group_by(Tratamento, ano) |>
  dplyr::summarise(
    media = mean(Rugosity_SD, na.rm = TRUE),
    se = sd(Rugosity_SD, na.rm = TRUE) / sqrt(dplyr::n()),
    .groups = "drop"
  ) |>
  tidyr::pivot_wider(
    names_from = ano,
    values_from = c(media, se),
    names_sep = "_"
  ) |>
  dplyr::mutate(
    variacao_media = round((media_2024 - media_2022) / 2, 2),
    SE_variacao = round(sqrt(se_2022^2 + se_2024^2) / 2, 2)) |>
  dplyr::mutate(values = paste(variacao_media, "pm", SE_variacao)) |>
  dplyr::select(Tratamento, values) |>
  tidyr::pivot_wider(names_from = Tratamento, values_from = values)
# A tibble: 1 × 5
  `Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%` `Planting 100%`
  <chr>        <chr>        <chr>        <chr>        <chr>          
1 0.34 pm 0.19 0.39 pm 0.12 0.34 pm 0.24 0.33 pm 0.13 0.11 pm 0.05   

d) LAI Mean

2022

show code
LAI2022 <- lidar_wide |>
  dplyr::filter(ano==2022) |>
  with(lme4::lmer(LAI_mean ~ Tratamento + (1|Bloco) ) )
par(mfrow=c(1,2))
DHARMa::plotResiduals(LAI2022)
DHARMa::plotQQunif(LAI2022)

2024

show code
LAI2024 <- lidar_wide |>
  dplyr::filter(ano==2024) |>
  with(lme4::lmer(LAI_mean ~ Tratamento + (1|Bloco) ) )
par(mfrow=c(1,2))
DHARMa::plotResiduals(LAI2024)
DHARMa::plotQQunif(LAI2024)

Full model

show code
modeloLAI <- lidar_wide |>
  with(lme4::lmer(LAI_mean ~ Tratamento * ano + (1 | Bloco)))
par(mfrow=c(1,2))
DHARMa::plotResiduals(modeloLAI)
DHARMa::plotQQunif(modeloLAI)

Pairwise comparison

show code
emmeans::emmeans(modeloLAI, pairwise ~ Tratamento | ano) |> multcomp::cld()
ano = 2022:
 Tratamento    emmean    SE   df lower.CL upper.CL .group
 Nuclei 25%     0.447 0.159 13.2    0.105    0.789  1    
 Strips 25%     0.591 0.159 13.2    0.249    0.933  1    
 Nuclei 50%     0.781 0.159 13.2    0.439    1.123  1    
 Strips 50%     0.828 0.159 13.2    0.485    1.170  1    
 Planting 100%  1.565 0.159 13.2    1.223    1.907   2   

ano = 2024:
 Tratamento    emmean    SE   df lower.CL upper.CL .group
 Nuclei 25%     0.900 0.159 13.2    0.558    1.242  1    
 Strips 25%     1.047 0.159 13.2    0.705    1.389  12   
 Strips 50%     1.353 0.159 13.2    1.011    1.695  12   
 Nuclei 50%     1.383 0.159 13.2    1.041    1.725   2   
 Planting 100%  2.408 0.159 13.2    2.066    2.750    3  

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 5 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Trend:

show code
emmeans::emtrends(modeloLAI, pairwise ~ Tratamento, var = "ano", poly.trend = TRUE)|> multcomp::cld()
 Tratamento    ano.trend     SE df lower.CL upper.CL .group
 Nuclei 25%        0.226 0.0812 36   0.0617    0.391  1    
 Strips 25%        0.228 0.0812 36   0.0633    0.393  1    
 Strips 50%        0.263 0.0812 36   0.0980    0.428  1    
 Nuclei 50%        0.301 0.0812 36   0.1360    0.466  1    
 Planting 100%     0.421 0.0812 36   0.2565    0.586  1    

Results are averaged over the levels of: ano 
Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 5 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Fixed vs random

show code
MuMIn::r.squaredGLMM(modeloLAI)
           R2m       R2c
[1,] 0.7064328 0.8459851
show code
0.7064328 /0.8459851
[1] 0.8350417
show code
sjPlot::tab_model(modeloLAI)
  Dependent variable
Predictors Estimates CI p
(Intercept) -457.36 -790.02 – -124.71 0.008
Tratamento [Strips 25%] -3.08 -473.52 – 467.36 0.989
Tratamento [Nuclei 50%] -150.02 -620.46 – 320.42 0.522
Tratamento [Strips 50%] -73.12 -543.56 – 397.32 0.755
Tratamento [Planting
100%]
-392.88 -863.32 – 77.56 0.099
ano 0.23 0.06 – 0.39 0.008
Tratamento [Strips 25%] ×
ano
0.00 -0.23 – 0.23 0.989
Tratamento [Nuclei 50%] ×
ano
0.07 -0.16 – 0.31 0.521
Tratamento [Strips 50%] ×
ano
0.04 -0.20 – 0.27 0.753
Tratamento [Planting
100%] × ano
0.19 -0.04 – 0.43 0.098
Random Effects
σ2 0.07
τ00 Bloco 0.06
ICC 0.48
N Bloco 5
Observations 50
Marginal R2 / Conditional R2 0.706 / 0.846
show code
broom.mixed::tidy(modeloLAI) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))
# A tibble: 12 × 6
   effect   group    term                        estimate std.error statistic
   <chr>    <chr>    <chr>                          <dbl>     <dbl>     <dbl>
 1 fixed    <NA>     (Intercept)                 -457.      164.       -2.78 
 2 fixed    <NA>     TratamentoStrips 25%          -3.08    232.       -0.013
 3 fixed    <NA>     TratamentoNuclei 50%        -150.      232.       -0.646
 4 fixed    <NA>     TratamentoStrips 50%         -73.1     232.       -0.315
 5 fixed    <NA>     TratamentoPlanting 100%     -393.      232.       -1.69 
 6 fixed    <NA>     ano                            0.226     0.081     2.79 
 7 fixed    <NA>     TratamentoStrips 25%:ano       0.002     0.115     0.014
 8 fixed    <NA>     TratamentoNuclei 50%:ano       0.074     0.115     0.647
 9 fixed    <NA>     TratamentoStrips 50%:ano       0.036     0.115     0.316
10 fixed    <NA>     TratamentoPlanting 100%:ano    0.195     0.115     1.70 
11 ran_pars Bloco    sd__(Intercept)                0.245    NA        NA    
12 ran_pars Residual sd__Observation                0.257    NA        NA    

Table S2

show code
lidar_wide |>
  dplyr::group_by(ano, Tratamento)|>
  dplyr::summarise(Mean = round(mean(LAI_mean),2),
                   SE = round((sd(LAI_mean)/sqrt(dplyr::n())), 2), 
                   .groups = "drop") |>
  dplyr::mutate(values = paste(Mean, "pm", SE)) |>
  dplyr::select(-c(Mean, SE)) |>
  tidyr::pivot_wider(names_from = Tratamento, values_from = values)
# A tibble: 2 × 6
    ano `Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%` `Planting 100%`
  <dbl> <chr>        <chr>        <chr>        <chr>        <chr>          
1  2022 0.45 pm 0.09 0.59 pm 0.12 0.78 pm 0.08 0.83 pm 0.16 1.57 pm 0.17   
2  2024 0.9 pm 0.15  1.05 pm 0.2  1.38 pm 0.16 1.35 pm 0.19 2.41 pm 0.21   

Increment/year (slope)

show code
lidar_wide |>
  dplyr::group_by(Tratamento, ano) |>
  dplyr::summarise(
    media = mean(LAI_mean, na.rm = TRUE),
    se = sd(LAI_mean, na.rm = TRUE) / sqrt(dplyr::n()),
    .groups = "drop"
  ) |>
  tidyr::pivot_wider(
    names_from = ano,
    values_from = c(media, se),
    names_sep = "_"
  ) |>
  dplyr::mutate(
    variacao_media = round((media_2024 - media_2022) / 2, 2),
    SE_variacao = round(sqrt(se_2022^2 + se_2024^2) / 2, 2)) |>
  dplyr::mutate(values = paste(variacao_media, "pm", SE_variacao)) |>
  dplyr::select(Tratamento, values) |>
  tidyr::pivot_wider(names_from = Tratamento, values_from = values)
# A tibble: 1 × 5
  `Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%` `Planting 100%`
  <chr>        <chr>        <chr>        <chr>        <chr>          
1 0.23 pm 0.09 0.23 pm 0.12 0.3 pm 0.09  0.26 pm 0.12 0.42 pm 0.14   

e) Figure - LiDAR

(i) Canopy height model (CHM)

show code
max_year <- lidar_wide |>
  dplyr::group_by(ano, Tratamento)|>
  dplyr::summarise(Mean = mean(CHM_mean),.groups="drop") |>
  dplyr::group_by(ano)|>
  dplyr::summarise(max_year = max(Mean),.groups="drop")
lidar_wide |>
  dplyr::group_by(ano, Tratamento)|>
  dplyr::summarise(Mean = mean(CHM_mean),.groups="drop") |>
  dplyr::right_join(max_year, by="ano") |>
  dplyr::mutate(porcentos = round( (Mean/max_year)*100, 0))
# A tibble: 10 × 5
     ano Tratamento     Mean max_year porcentos
   <dbl> <fct>         <dbl>    <dbl>     <dbl>
 1  2022 Nuclei 25%     1.89     4.91        39
 2  2022 Strips 25%     2.12     4.91        43
 3  2022 Nuclei 50%     3.00     4.91        61
 4  2022 Strips 50%     2.93     4.91        60
 5  2022 Planting 100%  4.91     4.91       100
 6  2024 Nuclei 25%     2.76     6.99        39
 7  2024 Strips 25%     3.16     6.99        45
 8  2024 Nuclei 50%     4.70     6.99        67
 9  2024 Strips 50%     4.43     6.99        63
10  2024 Planting 100%  6.99     6.99       100

Plot

show code
letras <- data.frame(
  Tratamento = rep(levels(as.factor(lidar_wide$Tratamento)),each=2),
  ano = c(2022, 2024),
  label = c("A      ","       a",
            "A      ","       a", 
            "A      ","       b",
            "A      ","       b",
            "B      ","       c"),
  y = c(2.8, 3.9, 
        3.1, 4.5, 
        3.9, 5.7, 
        3.9, 5.7, 
        5.9, 7.8)  )
porcentos <- data.frame(
  Tratamento = rep(levels(as.factor(lidar_wide$Tratamento)),each=2),
  ano = c(2022, 2024),
  label = c("39%        ","          39%",
            "43%        ","          45%",
            "61%        ","          67%",
            "60%        ","          63%",
            "100%         ","          100%"),  
  y = c(2.2, 3.1, 
        2.5, 3.7, 
        3.3, 5, 
        3.3, 5, 
        5.3, 7.1)            
)
p1 <- lidar_wide |>
  dplyr::group_by(ano, Tratamento)|>
  dplyr::summarise(Mean = mean(CHM_mean),
                   SD = sd(CHM_mean),
                   SE = SD/sqrt(dplyr::n()),
                   .groups="drop") |>
  ggplot(aes(x = Tratamento, y = Mean, fill = as.factor(ano))) +
  geom_bar(position = position_dodge2(preserve = "single"), 
           stat = "identity", color = "black") +
  geom_errorbar(aes(ymin = Mean - SE, ymax = Mean + SE), width = .2, position = position_dodge(0.9)) +
  scale_fill_manual(
    labels = c('5 years', '7 years'), 
    values = c("#CAF0F8", "#00B4D8")) +
  geom_text(data = letras, aes(x = Tratamento, y = y, label = label),size = 3.5, color="black", vjust = -0.5) +
  geom_text(data = porcentos, aes(x = Tratamento, y = y, label = label),size = 3, color="black", vjust = -0.5) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 9), breaks = seq(0, 8, by = 1))+
  labs(fill = "", 
       y =  "Canopy height model (m)",
       x = "") +
  theme_classic()+
    theme(
      axis.text.x = element_blank(),
      axis.text.y = element_text(color="black", size=9),
      axis.title.y = element_text(color="black", size=10),
      legend.position = c(0.15, 0.9),
      legend.text = element_text(size = 10),
      legend.key.size = unit(3, "mm"),
      plot.margin = margin(t=5, 0, 0, 0, "mm"))
p1

(ii) Canopy cover

show code
max_year <- lidar_wide |>
  dplyr::group_by(ano, Tratamento)|>
  dplyr::summarise(Mean = mean(Canopy_cover),.groups="drop") |>
  dplyr::group_by(ano)|>
  dplyr::summarise(max_year = max(Mean),.groups="drop")
lidar_wide |>
  dplyr::group_by(ano, Tratamento)|>
  dplyr::summarise(Mean = mean(Canopy_cover),.groups="drop") |>
  dplyr::right_join(max_year, by="ano") |>
  dplyr::mutate(porcentos = round( (Mean/max_year)*100, 0))
# A tibble: 10 × 5
     ano Tratamento     Mean max_year porcentos
   <dbl> <fct>         <dbl>    <dbl>     <dbl>
 1  2022 Nuclei 25%     31.6     79.1        40
 2  2022 Strips 25%     23.7     79.1        30
 3  2022 Nuclei 50%     26.5     79.1        33
 4  2022 Strips 50%     31.3     79.1        40
 5  2022 Planting 100%  79.1     79.1       100
 6  2024 Nuclei 25%     56.5     96.2        59
 7  2024 Strips 25%     44.7     96.2        46
 8  2024 Nuclei 50%     70.9     96.2        74
 9  2024 Strips 50%     69.8     96.2        73
10  2024 Planting 100%  96.2     96.2       100

plot

show code
letras <- data.frame(
  Tratamento = rep(levels(as.factor(lidar_wide$Tratamento)),each=2),
  ano = c(2022, 2024),
  label = c("A     ","       a",
            "A     ","       a", 
            "A     ","      ab",
            "A     ","      ab",
            "B     ","       b"),  
  y = c(52, 77, 
        40, 65, 
        40, 85, 
        45, 90, 
        97, 109)            
)
porcentos <- data.frame(
  Tratamento = rep(levels(as.factor(lidar_wide$Tratamento)),each=2),
  ano = c(2022, 2024),
  label = c("40%        ","          59%",
            "30%        ","          46%",
            "33%        ", "         74%",
            "40%        ","          73%",
            "100%         ","          100%"),  
  y = c(42, 67, 
        30, 55, 
        30, 75, 
        35, 80, 
        87, 100)  
)
p2 <- lidar_wide |>
  dplyr::group_by(ano, Tratamento)|>
  dplyr::summarise(Mean = mean(Canopy_cover),
                   SD = sd(Canopy_cover),
                   SE = SD/sqrt(dplyr::n()),
                   .groups="drop") |>
  ggplot(aes(x = Tratamento, y = Mean, fill = as.factor(ano))) +
  geom_bar(position = position_dodge2(preserve = "single"), 
           stat = "identity", color = "black") +
  geom_errorbar(aes(ymin = Mean - SE, ymax = Mean + SE), width = .2, position = position_dodge(0.9)) +
  scale_fill_manual(
    labels = c('5 years', '7 years'), 
    values = c("#CAF0F8", "#00B4D8")) +
  geom_text(data = letras, aes(x = Tratamento, y = y, label = label),size = 3.5, color="black", vjust = -0.5) +
  geom_text(data = porcentos, aes(x = Tratamento, y = y, label = label),size = 3, color="black", vjust = -0.5) +
    scale_y_continuous(expand = c(0, 0), limits=c(0, 120), breaks = seq(0, 100, by = 20))+
  labs(fill = "", 
       y =  "Canopy cover (%)",
       x = "") +
  theme_classic()+
    theme(
    axis.text.x = element_blank(),
    axis.text.y = element_text(color="black", size=10),
    axis.title.y = element_text(color="black", size=10),
    plot.margin = margin(t=5, 0, 0, 0, "mm"),
    legend.position = 'none')
p2

(iii) LAI

show code
max_year <- lidar_wide |>
  dplyr::group_by(ano, Tratamento)|>
  dplyr::summarise(Mean = mean(LAI_mean),.groups="drop") |>
  dplyr::group_by(ano)|>
  dplyr::summarise(max_year = max(Mean),.groups="drop")
lidar_wide |>
  dplyr::group_by(ano, Tratamento)|>
  dplyr::summarise(Mean = mean(LAI_mean),.groups="drop") |>
  dplyr::right_join(max_year, by="ano") |>
  dplyr::mutate(porcentos = round( (Mean/max_year)*100, 0))
# A tibble: 10 × 5
     ano Tratamento     Mean max_year porcentos
   <dbl> <fct>         <dbl>    <dbl>     <dbl>
 1  2022 Nuclei 25%    0.447     1.57        29
 2  2022 Strips 25%    0.591     1.57        38
 3  2022 Nuclei 50%    0.781     1.57        50
 4  2022 Strips 50%    0.828     1.57        53
 5  2022 Planting 100% 1.57      1.57       100
 6  2024 Nuclei 25%    0.900     2.41        37
 7  2024 Strips 25%    1.05      2.41        43
 8  2024 Nuclei 50%    1.38      2.41        57
 9  2024 Strips 50%    1.35      2.41        56
10  2024 Planting 100% 2.41      2.41       100

Plot

show code
letras <- data.frame(
  Tratamento = rep(levels(as.factor(lidar_wide$Tratamento)),each=2),
  ano = c(2022, 2024),
  label = c("A       ","       a",
            "A       ","      ab",
            "A       ", "      b",
            "A       ","      ab",
            "B       ","       c"),  
  y = c( 0.85, 1.4, 
         1.05, 1.55, 
         1.15, 1.82, 
         1.25, 1.82, 
         2.05, 2.9)
  ) |>
  dplyr::mutate(
    Tratamento = forcats::fct_recode(Tratamento,
                                     "Nuclei\n25%" = "Nuclei 25%",
                                     "Strips\n25%" = "Strips 25%",
                                     "Nuclei\n50%" = "Nuclei 50%",
                                     "Strips\n50%" = "Strips 50%",
                                     "Planting\n100%" = "Planting 100%"
                                     )
  )

porcentos <- data.frame(
  Tratamento = rep(levels(as.factor(lidar_wide$Tratamento)),each=2),
  ano = c(2022, 2024),
  label = c("29%        ","          37%",
            "38%       ","          43%",
            "50%       ","          57%",
            "53%       ","          56%",
            "100%          ","        100%"),  
  y = c( 0.55, 1.1, 
         0.75, 1.25, 
         0.85, 1.52, 
         0.95, 1.52, 
         1.75, 2.65)
  )|>
  dplyr::mutate(
    Tratamento = forcats::fct_recode(Tratamento,
                                     "Nuclei\n25%" = "Nuclei 25%",
                                     "Strips\n25%" = "Strips 25%",
                                     "Nuclei\n50%" = "Nuclei 50%",
                                     "Strips\n50%" = "Strips 50%",
                                     "Planting\n100%" = "Planting 100%"
                                     )
  )

p3 <- lidar_wide |>
  dplyr::group_by(ano, Tratamento)|>
  dplyr::summarise(Mean = mean(LAI_mean),
                   SD = sd(LAI_mean),
                   SE = SD/sqrt(dplyr::n()),
                   .groups="drop") |>
  dplyr::mutate(
    Tratamento = forcats::fct_recode(Tratamento,
                                     "Nuclei\n25%" = "Nuclei 25%",
                                     "Strips\n25%" = "Strips 25%",
                                     "Nuclei\n50%" = "Nuclei 50%",
                                     "Strips\n50%" = "Strips 50%",
                                     "Planting\n100%" = "Planting 100%"))|>
  ggplot(aes(x = Tratamento, y = Mean, fill = as.factor(ano))) +
  geom_bar(position = position_dodge2(preserve = "single"), 
           stat = "identity", color = "black") +
  geom_errorbar(aes(ymin = Mean - SE, ymax = Mean + SE), width = .2, position = position_dodge(0.9)) +
  scale_fill_manual(
    labels = c('5 years', '7 years'), 
    values = c("#CAF0F8", "#00B4D8")) +
  geom_text(data = letras, aes(x = Tratamento, y = y, label = label),size = 3.5, color="black", vjust = -0.5) +
  geom_text(data = porcentos, aes(x = Tratamento, y = y, label = label),size = 3, color="black", vjust = -0.5) +
  scale_y_continuous(expand = c(0, 0), limits = c(0,3.2), breaks = seq(0, 3, by=1))+
  labs(fill = "", 
       y =  expression("LAI (m"^{2} ~ "/ m"^{2}~ ")"),
       x = "") +
  theme_classic()+
  theme(
    axis.text.x = element_text(color="black", hjust = 0.5, size=10),
    axis.text.y = element_text(color="black", size=10),
    axis.title.y = element_text(color="black", size=10),
    plot.margin = margin(t=5, 0, 0, 0, "mm"),
    legend.position = 'none')

p3

(iv) All together

show code
# 90 mm x 180 mm
fig_lidar <- cowplot::plot_grid(p1,p2,p3, ncol = 1,align = "v", axis = "l",
       labels = c("(a)", "(b)", "(c)"), label_size = 10, label_x = 0.01, label_y = 1)

fig_lidar

show code
# ggsave("fig_lidar.png", plot=fig_lidar, width=90, height = 210 ,unit="mm")

capital letters: 5 years

small letters: 7 years

III. Results - Costs

A) Cost per unit (Biomass and Canopy coverage)

Table of Costs (in R$)

show code
costs <- data.frame(
  tratamento = c("Planting 100%",   "Strips 50%",   "Nuclei 50%",   "Strips 25%",   "Nuclei 25%"),
  labor_machinery =  c(7860.00, 6407.90, 7563.90, 5131.05, 6395.23),
  inputs = c(7032.50, 4096.25, 4096.25, 2628.13, 2628.13)) |>
  dplyr::mutate(total = labor_machinery + inputs)

Using raw values

show code
cover_raw <- lidar_wide |> 
  dplyr::select(ano, Bloco, Tratamento, Canopy_cover) |>
  dplyr::rename(bloco = Bloco, tratamento = Tratamento) |>
  dplyr::mutate(tratamento = ordered(tratamento, 
                                     levels = c("Nuclei 25%", "Strips 25%", "Nuclei 50%", "Strips 50%", "Planting 100%")),
                bloco = ordered(bloco, levels = c("BI","BII", "BIII", "BIV", "BV")),
                bloco = as.character(as.numeric(bloco)))|>
  dplyr::left_join(costs) |>
  dplyr::mutate(totalUS = (total/5)/5.8)
Joining with `by = join_by(tratamento)`
show code
cost_raw <- AGB_data |> 
  dplyr::select(ano, bloco, tratamento, AGB_Mgha) |> 
  dplyr::mutate(
    tratamento = dplyr::case_when(
      tratamento == "100" ~ "Planting 100%",
      tratamento == "F25" ~ "Strips 25%",
      tratamento == "F50" ~ "Strips 50%",
      tratamento == "N25" ~ "Nuclei 25%",
      tratamento == "N50" ~ "Nuclei 50%"
    )
    ) |>
  dplyr::right_join(cover_raw) |>
  dplyr::mutate(agb_cost = totalUS/AGB_Mgha,
                cc_cost = totalUS/Canopy_cover,
                ano = ano-2017) |>
  dplyr::select(ano:tratamento, agb_cost, cc_cost)
Joining with `by = join_by(ano, bloco, tratamento)`

(i) Tests

show code
modeloCC_cost <- lme4::lmer(log(cc_cost) ~ tratamento*ano + (1|bloco), data=cost_raw) 
par(mfrow=c(1,2))
DHARMa::plotResiduals(modeloCC_cost)
DHARMa::plotQQunif(modeloCC_cost)

show code
emmeans::emmeans(modeloCC_cost, pairwise ~ tratamento | ano) |> multcomp::cld()
ano = 5:
 tratamento    emmean    SE df lower.CL upper.CL .group
 Planting 100%   1.91 0.256 18     1.37     2.45  1    
 Strips 25%      2.50 0.256 18     1.96     3.03  12   
 Strips 50%      2.56 0.256 18     2.02     3.10  12   
 Nuclei 50%      2.73 0.256 18     2.19     3.27  12   
 Nuclei 25%      2.77 0.256 18     2.24     3.31   2   

ano = 7:
 tratamento    emmean    SE df lower.CL upper.CL .group
 Planting 100%   1.68 0.256 18     1.14     2.21  1    
 Strips 50%      1.75 0.256 18     1.21     2.29  1    
 Nuclei 50%      1.79 0.256 18     1.25     2.32  1    
 Strips 25%      1.92 0.256 18     1.38     2.46  1    
 Nuclei 25%      1.99 0.256 18     1.45     2.53  1    

Degrees-of-freedom method: kenward-roger 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 5 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
show code
sjPlot::tab_model(modeloCC_cost)
  log(cc cost)
Predictors Estimates CI p
(Intercept) 4.73 2.93 – 6.54 <0.001
tratamento [Nuclei 50%] 0.36 -2.15 – 2.86 0.775
tratamento [Planting
100%]
-2.23 -4.74 – 0.27 0.079
tratamento [Strips 25%] -0.81 -3.31 – 1.70 0.519
tratamento [Strips 50%] -0.15 -2.66 – 2.36 0.905
ano -0.39 -0.68 – -0.10 0.010
tratamento [Nuclei 50%] ×
ano
-0.08 -0.49 – 0.33 0.697
tratamento [Planting
100%] × ano
0.27 -0.14 – 0.69 0.186
tratamento [Strips 25%] ×
ano
0.11 -0.31 – 0.52 0.607
tratamento [Strips 50%] ×
ano
-0.01 -0.42 – 0.40 0.952
Random Effects
σ2 0.21
τ00 bloco 0.12
ICC 0.37
N bloco 5
Observations 50
Marginal R2 / Conditional R2 0.341 / 0.583
show code
broom.mixed::tidy(modeloCC_cost) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), round, 3)) 
# A tibble: 12 × 6
   effect   group    term                        estimate std.error statistic
   <chr>    <chr>    <chr>                          <dbl>     <dbl>     <dbl>
 1 fixed    <NA>     (Intercept)                    4.73      0.89      5.32 
 2 fixed    <NA>     tratamentoNuclei 50%           0.356     1.24      0.287
 3 fixed    <NA>     tratamentoPlanting 100%       -2.23      1.24     -1.80 
 4 fixed    <NA>     tratamentoStrips 25%          -0.806     1.24     -0.651
 5 fixed    <NA>     tratamentoStrips 50%          -0.149     1.24     -0.12 
 6 fixed    <NA>     ano                           -0.392     0.144    -2.72 
 7 fixed    <NA>     tratamentoNuclei 50%:ano      -0.08      0.204    -0.392
 8 fixed    <NA>     tratamentoPlanting 100%:ano    0.274     0.204     1.35 
 9 fixed    <NA>     tratamentoStrips 25%:ano       0.106     0.204     0.518
10 fixed    <NA>     tratamentoStrips 50%:ano      -0.012     0.204    -0.061
11 ran_pars bloco    sd__(Intercept)                0.347    NA        NA    
12 ran_pars Residual sd__Observation                0.455    NA        NA    
show code
modeloAGB_cost <- lme4::lmer(log(agb_cost) ~ tratamento*ano + (1|bloco), data=cost_raw ) 
par(mfrow=c(1,2))
DHARMa::plotResiduals(modeloAGB_cost)
DHARMa::plotQQunif(modeloAGB_cost)

show code
emmeans::emmeans(modeloAGB_cost, pairwise ~ tratamento | ano) |> multcomp::cld()
ano = 5:
 tratamento    emmean    SE   df lower.CL upper.CL .group
 Planting 100%   3.21 0.217 17.5     2.75     3.67  1    
 Strips 50%      3.68 0.217 17.5     3.22     4.14  1    
 Strips 25%      3.76 0.217 17.5     3.30     4.22  1    
 Nuclei 50%      3.81 0.217 17.5     3.36     4.27  12   
 Nuclei 25%      4.50 0.217 17.5     4.04     4.95   2   

ano = 7:
 tratamento    emmean    SE   df lower.CL upper.CL .group
 Planting 100%   2.79 0.217 17.5     2.33     3.25  1    
 Strips 25%      3.08 0.217 17.5     2.63     3.54  1    
 Strips 50%      3.12 0.217 17.5     2.67     3.58  12   
 Nuclei 50%      3.17 0.217 17.5     2.71     3.62  12   
 Nuclei 25%      3.80 0.217 17.5     3.35     4.26   2   

Degrees-of-freedom method: kenward-roger 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 5 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

And the slopes?

show code
emmeans::emtrends(modeloCC_cost, pairwise ~ tratamento, var = "ano", poly.trend = TRUE)|> multcomp::cld()
 tratamento    ano.trend    SE df lower.CL upper.CL .group
 Nuclei 50%       -0.472 0.144 36   -0.764 -0.18000  1    
 Strips 50%       -0.405 0.144 36   -0.697 -0.11250  1    
 Nuclei 25%       -0.392 0.144 36   -0.684 -0.10005  1    
 Strips 25%       -0.287 0.144 36   -0.579  0.00551  1    
 Planting 100%    -0.118 0.144 36   -0.410  0.17419  1    

Results are averaged over the levels of: ano 
Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 5 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
show code
emmeans::emtrends(modeloAGB_cost, pairwise ~ tratamento, var = "ano", poly.trend = TRUE)|> multcomp::cld()
 tratamento    ano.trend    SE df lower.CL upper.CL .group
 Nuclei 25%       -0.346 0.121 36   -0.592  -0.1006  1    
 Strips 25%       -0.338 0.121 36   -0.584  -0.0928  1    
 Nuclei 50%       -0.324 0.121 36   -0.570  -0.0782  1    
 Strips 50%       -0.280 0.121 36   -0.525  -0.0340  1    
 Planting 100%    -0.210 0.121 36   -0.456   0.0355  1    

Results are averaged over the levels of: ano 
Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 5 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
show code
sjPlot::tab_model(modeloAGB_cost)
  log(agb cost)
Predictors Estimates CI p
(Intercept) 6.23 4.71 – 7.74 <0.001
tratamento [Nuclei 50%] -0.79 -2.90 – 1.32 0.452
tratamento [Planting
100%]
-1.96 -4.07 – 0.15 0.067
tratamento [Strips 25%] -0.77 -2.88 – 1.34 0.463
tratamento [Strips 50%] -1.15 -3.26 – 0.96 0.278
ano -0.35 -0.59 – -0.10 0.007
tratamento [Nuclei 50%] ×
ano
0.02 -0.32 – 0.37 0.897
tratamento [Planting
100%] × ano
0.14 -0.21 – 0.48 0.432
tratamento [Strips 25%] ×
ano
0.01 -0.34 – 0.35 0.964
tratamento [Strips 50%] ×
ano
0.07 -0.28 – 0.41 0.700
Random Effects
σ2 0.15
τ00 bloco 0.09
ICC 0.38
N bloco 5
Observations 50
Marginal R2 / Conditional R2 0.498 / 0.688
show code
broom.mixed::tidy(modeloAGB_cost) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), round, 3)) 
# A tibble: 12 × 6
   effect   group    term                        estimate std.error statistic
   <chr>    <chr>    <chr>                          <dbl>     <dbl>     <dbl>
 1 fixed    <NA>     (Intercept)                    6.23      0.749     8.31 
 2 fixed    <NA>     tratamentoNuclei 50%          -0.792     1.04     -0.76 
 3 fixed    <NA>     tratamentoPlanting 100%       -1.96      1.04     -1.88 
 4 fixed    <NA>     tratamentoStrips 25%          -0.773     1.04     -0.742
 5 fixed    <NA>     tratamentoStrips 50%          -1.15      1.04     -1.1  
 6 fixed    <NA>     ano                           -0.346     0.121    -2.86 
 7 fixed    <NA>     tratamentoNuclei 50%:ano       0.022     0.171     0.13 
 8 fixed    <NA>     tratamentoPlanting 100%:ano    0.136     0.171     0.794
 9 fixed    <NA>     tratamentoStrips 25%:ano       0.008     0.171     0.045
10 fixed    <NA>     tratamentoStrips 50%:ano       0.067     0.171     0.389
11 ran_pars bloco    sd__(Intercept)                0.298    NA        NA    
12 ran_pars Residual sd__Observation                0.383    NA        NA    

(ii) Figure Cost-effectiveness

show code
fig_data_cost_raw <- cost_raw |>
  dplyr::group_by(ano, tratamento)|>
  dplyr::summarise(
    cost_cc = mean(cc_cost),
    SE_cc = sd(cc_cost)/sqrt(dplyr::n()),
    cost_agb = mean(agb_cost),
    SE_agb = sd(agb_cost)/sqrt(dplyr::n()), 
    .groups="drop")|>
  dplyr::mutate(tratamento = ordered(tratamento, levels = c("Nuclei 25%", "Strips 25%", "Nuclei 50%", "Strips 50%", "Planting 100%")),
                tratamento = forcats::fct_recode(tratamento,
                                     "Nuclei\n25%" = "Nuclei 25%",
                                     "Strips\n25%" = "Strips 25%",
                                     "Nuclei\n50%" = "Nuclei 50%",
                                     "Strips\n50%" = "Strips 50%",
                                     "Planting\n100%" = "Planting 100%"
                                     ))
show code
letras <- data.frame(
  tratamento = rep(levels(as.factor(fig_data_cost_raw$tratamento)), each=2),
  ano = c(5, 7),
  label = c("B      ","        a",
            "AB      ","       a", 
            "AB      ","       a",
            "AB      ","       a",
            "A      ","       a"),
  y = c(33, 16, 
        18, 10, 
        18, 8, 
        18, 10, 
        10, 6))

fig_cost1 <- fig_data_cost_raw |>
  ggplot(aes(x = tratamento, y = cost_cc, fill = as.factor(ano))) +
  geom_bar(position = position_dodge2(preserve = "single"), 
           stat = "identity", color = "black") +
  geom_errorbar(aes(ymin = cost_cc - SE_cc, ymax = cost_cc + SE_cc), width = .2, position = position_dodge(0.9)) +
  scale_fill_manual(
    labels = c('5 years', '7 years'), 
    values = c("#CAF0F8", "#00B4D8")) +
  geom_text(data = letras, aes(x = tratamento, y = y, label = label), size = 3.5, color="black", vjust = -0.5) +
  scale_y_continuous(expand=c(0,0), limits=c(0,38), breaks = c(0, 10,20, 30))+
  labs(fill = "",
       y = expression("CE Canopy cover (" * "\u0024" * ".%"^{-1} * ")"),
       x = "") +
  theme_classic()+
    theme(
    axis.text.x = element_blank(),
    axis.text.y = element_text(color="black", size=10),
    axis.title.y = element_text(color="black", size=10),
    plot.margin = margin(t=5, 0, 0, 0, "mm"),
    legend.position = c(0.85, 0.95))

fig_cost1

show code
letras <- data.frame(
  tratamento = rep(levels(as.factor(fig_data_cost_raw$tratamento)), each=2),
  ano = c(5, 7),
  label = c("B     ","       b",
            "A     ","      a", 
            "AB    ","     ab",
            "A     ","     ab",
            "A     ","      a"),
  y = c(210, 100, 53, 25, 53, 30, 53, 33, 30, 20)            
)


fig_cost2 <- fig_data_cost_raw |>
  ggplot(aes(x = tratamento, y = cost_agb, fill = as.factor(ano))) +
  geom_bar(position = position_dodge2(preserve = "single"), 
           stat = "identity", color = "black") +
  geom_errorbar(aes(ymin = cost_agb - SE_agb, ymax = cost_agb + SE_agb), 
                width = .2, position = position_dodge(0.9)) +
  scale_fill_manual(
    labels = c('5 years', '7 years'), 
    values = c("#CAF0F8", "#00B4D8")) +
  geom_text(data = letras, aes(x = tratamento, y = y, label = label), size = 3.5, color="black", vjust = -0.5) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 240), breaks = c(0,50,100, 150, 200))+
  labs(fill = "",
       y = expression("CE AGB (" * "\u0024" * ".Mg"^{-1} * ")"),
       x = "") +
  theme_classic()+
    theme(
    axis.line.x = element_line(color = "black", linewidth = 0.8, linetype = "solid"),
    axis.line.y = element_line(color = "black", linewidth = 0.8, linetype = "solid"),
    axis.ticks = element_line(color = "black", linewidth = 0.5),
    axis.text.x = element_text(color="black", hjust = .5, size=10),
    axis.text.y = element_text(color="black", size=10),
    legend.position = 'none',
    axis.title.y = element_text(color="black", size=10))
fig_cost2

show code
fig_cost <- cowplot::plot_grid(fig_cost1, fig_cost2, ncol = 1,
                               align = "v", axis = "l",
                               labels = c("(a)", "(b)"), label_size = 8, 
                               label_x = 0.05, label_y = c(1,1.1))
fig_cost

show code
# ggsave("fig_cost.png", plot=fig_cost, width=90, height = 120 ,unit="mm")

(iii) Table Cost-effectiveness

show code
fig_data_cost_raw |>
  dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 2))) 
# A tibble: 10 × 6
     ano tratamento       cost_cc SE_cc cost_agb SE_agb
   <dbl> <ord>              <dbl> <dbl>    <dbl>  <dbl>
 1     5 "Nuclei\n25%"      23.3   9.43    133.   69.9 
 2     5 "Nuclei\n50%"      15.5   1.08     46.4   4.84
 3     5 "Planting\n100%"    7.07  1.09     25.6   3.28
 4     5 "Strips\n25%"      12.8   1.8      44.6   5.95
 5     5 "Strips\n50%"      14.5   3.44     43.3   8.74
 6     7 "Nuclei\n25%"      10.5   4.82     58.6  26.0 
 7     7 "Nuclei\n50%"       6.26  0.94     24.4   2.88
 8     7 "Planting\n100%"    5.35  0.11     16.6   1.47
 9     7 "Strips\n25%"       7.52  1.37     22.6   2.66
10     7 "Strips\n50%"       6.63  1.92     26.3   7.65